Investigating the Origin of Azimuthal Site Response Variability using KiK-net Data: A Twofold Spectral Ratio Approach

Investigating the Origin of Azimuthal Site Response Variability using KiK-net Data: A Twofold Spectral Ratio Approach | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Research Article Investigating the Origin of Azimuthal Site Response Variability using KiK-net Data: A Twofold Spectral Ratio Approach Masashi Aoki, Yu Yamamoto, Ryoichi Tokumitsu, Yasuo Uchiyama, and 1 more This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-7806012/v1 This work is licensed under a CC BY 4.0 License Status: Posted Version 1 posted You are reading this latest preprint version Abstract In this study, we show that the widely observed azimuthal dependence of the surface-to-borehole spectral ratios (SBSRs) is governed primarily by local, site-specific complexities rather than regional path-specific effects, resolving a fundamental ambiguity in site response analysis. This ambiguity arises from the significant event-to-event variability seen in empirical data, which contradicts the stable, time-invariant theoretical transfer function predicted by one-dimensional (1D) site response analysis and introduces considerable epistemic uncertainty into ground motion models (GMMs). A key challenge has been the indeterminate origin of this variance whether attributable to the propagation path or the local site. To systematically decompose these contributing effects, we developed an analytical framework based on the twofold spectral ratio (TSR) method. Applying this framework to a large dataset from 293 nearby vertical array station pairs from Japan's Kiban–Kyoshin network (KiK-net), our analysis was designed to isolate and nullify the common path-effect term, enabling a direct, quantitative evaluation of site-specific variability. The analysis yielded two empirical results. First, the path-effect term isolated by the TSR showed high correlation between surface and borehole records, with the median correlation coefficient remaining above 0.5 across the entire 1–10 Hz frequency band, validating the successful separation of this component. Second, the residual event-to-event SBSR fluctuations between adjacent stations exhibited a near-zero spatial correlation across the entire frequency range, which is inconsistent with the hypothesis of a shared, regional cause. This reclassification of variability from a random uncertainty to a deterministic, site-specific attribute has direct implications for advancing seismic hazard assessment. Our dual-dataset approach provides a framework to validate the decomposition of site-specific variability into its aleatory and epistemic components. This provides the physical basis for developing non-ergodic GMMs, which can reduce overall uncertainty by explicitly modeling the predictable epistemic component, facilitating more sophisticated, physics-based site characterization. Surface-to-Borehole Spectral Ratio Azimuthal Variation Twofold Spectral Ratio Site Amplification Site-Specific Variability KiK-net Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 1. Introduction Site-specific ground motion amplification is a critical factor in earthquake-resistant design and seismic hazard assessment. One-dimensional (1D) site response analysis has long been a fundamental tool in geotechnical earthquake engineering practice, simplifying wave propagation to vertically propagating shear waves through horizontal soil layers (Kramer, 1996). Within this framework, site amplification is characterized by a theoretical transfer function (TTF). The TTF assumes ergodicity, where a site’s response is expected to be stable and time-invariant, representing a unique and predictable property of its local geology. A key metric for empirically measuring the TTF is the surface-to-borehole spectral ratio (SBSR). Pioneering studies demonstrated that by analyzing simultaneous recordings from surface and subsurface instruments, it was possible to compute actual amplification functions—often treated as the empirical transfer function (ETF) of the site—and validate the core principles of site response (Shima, 1962 ; Joyner et al., 1976 ). Extensive vertical array networks, such as the Kiban–Kyoshin network (KiK-net) in Japan, the California Strong-Motion Instrumentation Program (CSMIP), and the Strong-Motion Downhole Array in the Taipei Basin, have provided a substantial volume of data, facilitating the widespread application of the SBSR method. Nevertheless, the interpretation of the SBSR as a pure transfer function is complicated. Ideally, the SBSR would capture a true resonance of the site, which corresponds to the amplification at its natural frequencies under outcrop conditions. However, even within an idealized 1D framework, the interference from downward-propagating waves creates spectral troughs in the borehole recording, which can result in pseudoresonance peaks in the calculated SBSR. In practice, this theoretical complexity is further compounded, as observed borehole seismograms are also contaminated by complex wavefields, such as scattered waves from 3D geological irregularities (Thompson et al., 2012 ; Tao and Rathje, 2020 ). These physical complexities result in significant event-to-event variability, a phenomenon that directly challenges the ergodic assumption of 1D models, even for weak motions within the linear response range (Archuleta et al., 1992 ; Pilz and Cotton, 2019 ). The overall variability reflects the limitations of 1D models, contributing to epistemic uncertainty in ground motion models (GMMs) and probabilistic seismic hazard analysis (PSHA) (Parker et al., 2021 ). The event-to-event variability in site responses can be attributed to physical factors including the wave type, back-azimuth, and incidence angle (Field and Jacob, 1995 ). The presence of azimuth-dependent variability is a critical issue for modern data-driven methods that aim to break the ergodic assumption by separating local site amplification from regional path effects. This gives rise to a fundamental ambiguity, as the physical origin of this azimuthal dependence—whether it is attributable to the site term or the path term—has remained unresolved. Quantitatively distinguishing between these two sources is a crucial step for advancing site characterization. The practical importance of this distinction is evident in the development of non-ergodic GMMs. For instance, Ikeura ( 2020 ) proposed a method to investigate path characteristics without assuming a uniform attenuation function by first removing site effects, which were quantified as relative site factors (RSFs). The stability and validity of such an RSF, however, is directly challenged if it exhibits an azimuthal dependence of unknown origin. This highlights the need to resolve whether the dependence is a site-related or path-related phenomenon. The objective of this study is to investigate the physical origin of the observed azimuthal dependence of SBSRs. We employ an analytical approach based on the twofold spectral ratio (TSR) method, adapted for pairs of nearby vertical array stations, to systematically decompose site response variability into path-related and site-related components. 2. Data and Preliminary Analysis The dataset and preprocessing steps are described in this section. It also presents a preliminary analysis to characterize SBSR variability and to identify the primary factors influencing this variability. 2.1. KiK-net Strong-Motion Data Data were obtained from the KiK-net, a nationwide strong-motion seismograph network in Japan operated by the National Research Institute for Earth Science and Disaster Resilience (NIED). The network comprises approximately 700 stations, each with surface and borehole three-component accelerometers. This configuration is ideal for directly observing site amplification effects and provides an extensive dataset for our research (NIED 2019). 2.2. Earthquake Record Selection To compile a comprehensive, high-quality dataset, we selected earthquake records on the basis of several criteria, with the primary objective of analyzing the linear or quasilinear behavior of site responses. The selection criteria are as follows: Earthquake Magnitude : To ensure a sufficient signal-to-noise ratio for stable spectral analysis, only earthquakes with a Japan Meteorological Agency (JMA) -defined magnitude (M J ) of 4.0 or greater were included. Peak Ground Acceleration (PGA) : Records were selected if the surface PGA varied between 1.0 and 50.0 cm/s². The upper limit was intended to exclude records that may have been significantly affected by soil nonlinearity, following the approach of Noguchi and Sasatani ( 2011 ). Event Uniqueness : To isolate the ground motion from a single source, records were screened to exclude instances in which multiple earthquakes occurred in close succession by examining theoretical travel-time predictions via the JMA2001 travel-time table (Ueno et al. 2002 ). After applying these criteria, our final dataset comprised 152,352 records from 6,973 earthquakes observed at 693 stations. The distribution of key dataset parameters is shown in Fig. 1, which includes histograms of M J , focal depth, surface PGA, epicentral distance, hypocentral distance, and geometric incidence angle. The geometric incidence angle is the angle calculated from the straight-line path between the hypocenter and the station, ignoring reflections and refractions at layer boundaries. 2.3. Data Preprocessing All acceleration time histories were subjected to a series of standard preprocessing steps. Baseline correction was performed by subtracting the mean of the entire record. A critical step in analyzing borehole records is correcting the orientation of horizontal components. In this study, we applied correction angle data published by Shiomi ( 2012 ). For the 33 stations lacking such data, an independent orientation correction was performed by maximizing the coherence between the surface and borehole horizontal recordings in a low-frequency band. To analyze the SH-wave response, the two horizontal components were rotated to the great-circle path from the epicenter to the station to obtain the transverse component. The SBSR was then computed by dividing the Fourier amplitude spectrum of the transverse-component record at the surface by that in the borehole. To obtain a stable spectral ratio and clarify the trend of site response, the resulting ratio was smoothed using a Parzen window (Noguchi and Sasatani, 2011 ). A bandwidth of 1.0 Hz was selected to enhance the stability of the spectral ratio across the frequency range of interest. 2.4. Preliminary Analysis of SBSR Variability To understand the nature of the observed site response, the preliminary analysis focused on the 1–10 Hz frequency range, which is a band of primary interest for many engineering applications (Kramer, 1996). The analysis results are shown in Fig. 2. The standard deviation (SD) of the natural logarithm of the SBSR (ln(SBSR)) for all 693 stations is shown in Fig. 2(a). The results for the station with the largest maximum SD (MYGH13) are highlighted in red. This panel reveals that the SBSR variability is generally small at low frequencies but increases significantly with frequency. The azimuthal dependence of the response is shown in Fig. 2(b), which shows the mean SBSR for station MYGH13, grouped into four back-azimuth quadrants. The spectral shape clearly differs for each direction. To objectively classify the observed SBSR peaks, a Gaussian mixture model—a probabilistic model used to represent subpopulations within an overall population—was applied using the scikit-learn library in Python (Pedregosa et al, 2011 ) (Fig. 2(c)). The number of clusters in the model was defined on the basis of the average number of peaks observed per event at a given station. In the subsequent statistical analysis, we focused on the cluster with the lowest mean frequency. The variability in the lowest-frequency cluster, which could be caused by either a true resonance or pseudoresonance peak resulting from downward-propagating waves (Thompson et al. 2012 ; Tao and Rathje 2020 ), was considered the fundamental feature in this analysis. The clustering results for a representative station are shown in Fig. 2(c). Finally, one-way multivariate analysis of variance (MANOVA) of the peaks belonging to the lowest-frequency cluster was conducted to quantitatively compare the influences of the back-azimuth and geometric incidence angle, implemented with the statsmodels library in Python (Seabold and Perktold, 2010 ). The results are shown in the scatter plot of Pillai’s trace statistics in Fig. 2(d). MANOVA is a statistical method for comparing the means of two or more populations and for determining whether these differences are statistically significant. This technique has also been applied in fields such as geostatistical modeling (Sgobba et al. 2024 ). Pillai's trace was selected as the test statistic for our MANOVA because it is particularly robust against violations of assumptions, such as unequal sample sizes between groups, making it suitable for analyzing earthquake datasets (Hair et al., 2018 ). Among the 693 total stations, 347 stations provided sufficient data for a valid analysis. The results revealed that the statistics for the back-azimuth are more notable than those for the geometric incidence angle at most stations. Specifically, the back-azimuth was the dominant factor at 303 stations (87.3%), whereas the geometric incidence angle dominated at only 44 (12.7%). This preliminary analysis empirically establishes that the back-azimuth is the primary factor controlling the observed SBSR variability, raising the critical question of whether this considerable azimuthal dependence is caused by complexities along the wave propagation path or by local site conditions. 3. Methodology: Modified Twofold Spectral Ratio Approach The preliminary analysis in the previous section established that the observed SBSR variability is controlled primarily by the back-azimuth of incoming seismic waves, necessitating a distinction between two potential sources for the observed azimuthal dependence, namely, path- and site-specific effects. The methodology presented herein is therefore designed to quantitatively separate these two effects for analysis. 3.1. Analytical Framework The TSR method was established to estimate path-dependent seismic wave characteristics, most notably the quality factor (Q value) for anelastic attenuation, by systematically canceling common source and site amplification effects from observed ground motion spectra (Matsuzawa et al. 2004 ; Izutani 2000). In this study, we adapted the TSR framework not to estimate path effects but to first isolate and then eliminate them, enabling a focused evaluation of site-specific contributions to SBSR variability. The Fourier spectrum of an observed ground motion record is typically modeled in the frequency domain as a convolution of source, path, and site effects. For a simplified point-source and far-field approximation, this can be expressed as a product. The observed Fourier spectra at the surface ( \(\:{F}_{m,i}^{S}\) ) and in the borehole ( \(\:{F}_{m,i}^{B}\) ) for station i and earthquake event m are thus defined as follows (Aki and Richards, 2002 ): $$\:{F}_{m,i}^{S}\left(f\right)={S}_{m}\left(f\right)\bullet\:{P}_{m,i}\left(f\right)\bullet\:{G}_{m,i}^{S}\left(f\right)$$ 1 $$\:{F}_{m,i}^{B}\left(f\right)={S}_{m}\left(f\right)\bullet\:{P}_{m,i}\left(f\right)\bullet\:{G}_{m,i}^{B}\left(f\right)$$ 2 where \(\:{S}_{m}\left(f\right)\) denotes the source spectrum for event m , \(\:{P}_{m,i}\left(f\right)\) denotes the path transfer function from source m to station site i , \(\:{G}_{m,i}^{S}\left(f\right)\) is the surface site amplification factor at station i , and \(\:{G}_{m,i}^{B}\left(f\right)\) is the corresponding site response term at the borehole depth. It is essential to recognize that the borehole recording does not constitute pure input motion. The wavefield at the borehole sensor is invariably contaminated by downward-propagating waves reflected from the free surface and any overlying low-velocity layers. This downgoing wave effect complicates the interpretation of borehole records for use as a reference for incident motion. The interference between upward- and downward-propagating waves can produce spectral troughs in borehole recordings, which in turn are manifested as apparent or pseudoresonance peaks in the calculated SBSR that are not representative of simple 1D resonance of the soil column (Thompson et al. 2012 ; Tao and Rathje 2020 ). The SBSR, which is the direct ETF between the surface and borehole, can be calculated as follows: $$\:{SBSR}_{m,i}\left(f\right)={F}_{m,i}^{S}\left(f\right)/{F}_{m,i}^{B}\left(f\right)={G}_{m,i}^{S}\left(f\right)/{G}_{m,i}^{B}\left(f\right)$$ 3 In this ratio, the shared source and path terms are theoretically eliminated, leaving a relatively stable characteristic of the vertical site structure, with a theoretical dependence on the incidence angle that is negligible for near-vertically propagating waves. The TSR method is applied to a pair of stations ( i and j ) and a pair of events ( m and n ). The TSR for surface recordings ( \(\:{TSR}^{S}\) ) is formulated by computing a ratio of spectral ratios. This operation cancels the source terms and, under the assumption of linear and event-invariant site responses, cancels the site terms, isolating the differential path effect: $$\:{TSR}^{S}\left(f\right)=\frac{{F}_{m,i}^{S}\left(f\right)/{F}_{n,i}^{S}\left(f\right)}{{F}_{m,j}^{S}\left(f\right)/{F}_{n,j}^{S}\left(f\right)}=\frac{{P}_{m,i}\left(f\right)/{P}_{n,i}\left(f\right)}{{P}_{m,j}\left(f\right)/{P}_{n,j}\left(f\right)}$$ 4 The same operation is applied to borehole records ( \(\:{TSR}^{B}\) ) to isolate the differential path-effect term: $$\:{TSR}^{B}\left(f\right)=\frac{{F}_{m,i}^{B}\left(f\right)/{F}_{n,i}^{B}\left(f\right)}{{F}_{m,j}^{B}\left(f\right)/{F}_{n,j}^{B}\left(f\right)}=\frac{{P}_{m,i}\left(f\right)/{P}_{n,i}\left(f\right)}{{P}_{m,j}\left(f\right)/{P}_{n,j}\left(f\right)}$$ 5 The central proposition of this study is that by computing the ratio of these two TSRs, i.e., the dual-TSR (DTR), common differential path effects, which both ratios are designed to capture, can be canceled. $$\:DTR\left(f\right)={TSR}^{S}\left(f\right)/{TSR}^{B}\left(f\right)$$ 6 As a direct result of this formulation, the resulting ratio is mathematically sensitive to any violation of the assumption that the empirical site response (i.e., the SBSR) is independent of the earthquake event. By substituting the definition of the SBSR in Eq. ( 3 ) into Eq. ( 6 ), the DTR can be rearranged into the following equivalent form: $$\:DTR\left(f\right)=\frac{{SBSR}_{m,i}\left(f\right)/{SBSR}_{n,i}\left(f\right)}{{SBSR}_{m,j}\left(f\right)/{SBSR}_{n,j}\left(f\right)}$$ 7 This final expression reveals the physical meaning of the DTR: it compares the interevent SBSR fluctuation at station i to that at station j . If these fluctuations are driven by common path effects, the fluctuations at both stations should be highly correlated, and \(\:DTR\left(f\right)\) should approach unity. Conversely, if the fluctuations are driven by site-specific factors (e.g., 2D/3D scattering effects dependent on the arrival direction), they will be largely independent between the two stations, and \(\:DTR\left(f\right)\) will deviate significantly from unity. 3.2. Analysis Procedure The selection of station pairs for the final analysis was governed by two criteria designed to balance statistical robustness with physical constraints. First, the horizontal interstation distance had to be less than 20 km. This threshold represents a practical balance between securing a robust number of pairs for statistical analysis and limiting significant decorrelation of the incoming wavefield. Second, to ensure statistical stability, we required that each station pair have a minimum of 20 common earthquake pairs for each of the two analysis datasets described below. After applying these filtering steps, a final analysis set of 293 station pairs was selected (Fig. 4(a)). To systematically investigate the cause of SBSR variability, two distinct datasets of earthquake pairs were established for TSR analysis. The criteria for these datasets were designed to test specific hypotheses regarding the origin of variability. An example of the geographical distribution of earthquakes used to obtain these pairs relative to the MYGH04–MYGH13 station pair is shown in Fig. 4(b) and 4(c). Collocated Earthquake Pairs (Dataset A) : This dataset comprises pairs of earthquakes that satisfy two conditions: a hypocenter separation smaller than 5 km and a M J difference smaller than 0.2. This strategy was designed to minimize variations in both the propagation path and source size, thereby isolating SBSR variability caused by differences in the source radiation pattern or the complex site response to slight variations in the incoming wavefield. An example of this earthquake selection is shown for the MYGH04–MYGH13 station pair in Fig. 4(b). Common-Distance, Different-Path Pairs (Dataset B) : This dataset includes pairs of earthquakes that satisfy three conditions: a difference in hypocentral distance to the station pair smaller than 5 km, a M J difference smaller than 0.2, and a back-azimuth separation of at least 45 degrees. This strategy was designed to compare events that travel over different paths but exhibit similar overall travel distances. These criteria for M J and hypocentral distance ensure that the source term and the geometric attenuation are comparable for both events in a pair. However, because the paths encompass different crustal volumes, the anelastic attenuation (Q value) is expected to differ. This dataset facilitates a direct assessment of whether the TSR method consistently captures these differing path effects at both the surface and the borehole, serving as validation of the overall approach. An example of this earthquake selection is shown for the MYGH04–MYGH13 station pair in Fig. 4(c). For each of the 293 station pairs and for each of the two earthquake pairing strategies, two correlation analyses were conducted across the 1–10 Hz frequency range. These analyses provided the quantitative metrics for this study. Path-Effect Correlation Pearson’s correlation coefficient between the natural logarithm of the surface-TSR (ln( \(\:{TSR}^{S}\) )) and the borehole-TSR (ln( \(\:{TSR}^{B}\) ​)) was calculated. This calculation serves as a validation or quality control step for the methodology. A high correlation coefficient indicates that the TSR method can effectively isolate the same path-related term at both the surface and the borehole. Site-Specific Variability Correlations Pearson’s correlation coefficient between the logarithmic interevent SBSR fluctuations at the two sites (between ln[ \(\:{SBSR}_{m,i}\left(f\right)/{SBSR}_{n,i}\left(f\right)\) ] and ln[ \(\:{SBSR}_{m,j}\left(f\right)/{SBSR}_{n,j}\left(f\right)\) ​]) was calculated. This represents the primary metric of the study. It directly reflects the degree of spatial coherence in the observed site response variability. A low correlation coefficient suggests that the event-to-event fluctuations in the SBSR at one station are independent of those at a nearby station, strongly suggesting that variability is controlled by local, site-specific conditions. 4. Results and Discussion This section presents the results of the modified TSR approach, quantitatively assessing the spatial coherence of event-to-event SBSR variability. The analysis is presented first for the overall statistical distributions of all 293 station pairs, followed by a detailed examination of a representative station pair (MYGH04–MYGH13) to illustrate the underlying physical mechanisms. 4.1. Path-Effect Correlation: Surface-TSR versus Borehole-TSR The correlation between the surface-TSR (ln( \(\:{TSR}^{S}\) )) and the borehole-TSR (ln( \(\:{TSR}^{B}\) )) serves as a critical validation metric for our approach. A high correlation indicates that the TSR successfully isolates the same differential path-effect term at both the surface and borehole sensors, regardless of the complex wavefield in the borehole. The frequency-dependent Pearson correlation coefficients for this path-effect term are shown in Fig. 5. The results for the collocated earthquake pairs (Dataset A) are provided in the lower panel of Fig. 5(a), while the results for the common-distance, different-path pairs (Dataset B) are shown in the lower panel of Fig. 5(b). For both datasets, the overall statistical distributions demonstrate strong positive correlations between the surface-TSR and the borehole-TSR. In the case of the collocated pairs in Fig. 5(a) (lower panel), the results of the analysis show that the median correlation coefficient is approximately 0.9 at 1 Hz, and while it gradually decreases with increasing frequency, it remains above 0.5 up to 10 Hz. This finding indicates that for more than half of the station pairs, the path-related spectral modulations are captured consistently at both sensor depths. The results for the different-path dataset in Fig. 5(b) (lower panel) reveal a very similar trend, with a robust median correlation across the entire frequency band. This high degree of correlation across a large and diverse set of station pairs confirms the fundamental premise of the analysis: the path-effect term, \(\:{P}_{m,i}\left(f\right)/{P}_{n,i}\left(f\right)\) , is largely independent of the local site effect and is consistently captured at both the surface and the borehole. This consistency holds even when the paths substantially differ, as in Dataset B, which suggests that the anelastic attenuation characteristics (related to the Q value) along different paths are robustly captured by the TSR. The successful isolation of this common path term is essential, as it allows us to attribute the remaining variability to other factors. 4.2. Site-Specific Variability Correlation: Interstation SBSR Fluctuations Having validated the method's ability to isolate a common path term in Section 4.1 , the core of our subsequent analysis lies in determining whether the event-to-event SBSR variability is a spatially coherent effect driven by regional path phenomena, or a spatially incoherent effect driven by local site conditions. We therefore calculated the correlation of the interevent SBSR fluctuations ( \(\:{SBSR}_{m,i}\left(f\right)/{SBSR}_{n,i}\left(f\right)\) ) between each of the 293 station pairs. The results are shown in the upper panels of Fig. 5(a) and 5(b) and are consistent across both datasets. The median interstation correlation of the SBSR fluctuations remains close to zero across the entire 1–10 Hz frequency range. The interquartile range (IQR), indicated by the dark shaded area, is also narrowly centered on zero, typically ranging from approximately − 0.3 to + 0.3. This near-zero correlation indicates a lack of spatial coherence in SBSR variability. If the azimuthal dependence of the SBSR was caused primarily by large-scale, path-dependent effects such as directional variations in regional attenuation or scattering from deep crustal heterogeneities, one would expect to observe a positive correlation between nearby stations, as they would share a significant portion of their propagation paths. The absence of such a correlation strongly suggests that the observed event-to-event fluctuations are controlled primarily by factors unique to each individual site. These site-specific complexities, such as small-scale geological irregularities, surface topography, or basin-edge effects, are known to interact with the incoming wavefield in a manner that is highly sensitive to the arrival direction (Bard and Bouchon, 1985 ). The resulting scattered wavefield and directional resonance patterns are therefore unique to the immediate vicinity of each station, leading to the observed spatial decorrelation of SBSR variability even between stations separated by less than 20 km. 4.3. Case Study: Influence of Azimuthal Dependence at MYGH04–MYGH13 The overall statistical results in Section 4.2 indicate that SBSR variability is a spatially incoherent phenomenon controlled by local site conditions. To explore the physical mechanisms underlying this general conclusion in greater detail, we examined the results for the MYGH04–MYGH13 station pair, highlighted by magenta lines in Fig. 5. As established in the preliminary analysis (Fig. 2(b)), the MYGH13 station exhibits exceptionally high azimuthal dependence of its SBSR. Let us first consider the results for Dataset B (Fig. 5(b)), which represents the more complex scenario. The path-effect correlation for this pair (lower panel, magenta line) deviates significantly from the median trend. It indicates notable decreases in correlation, decreasing to nearly zero within the 3–4.5 Hz and 6–9.5 Hz frequency bands. These frequency bands coincide precisely with the frequencies at which the SBSR at MYGH13 exhibits its largest peak shifts and shape changes with back-azimuth. This result is instructive: it demonstrates that when site effects highly depend on the arrival direction, they are no longer completely canceled during TSR calculation. The site term becomes event-dependent (i.e., \(\:{G}_{m,i}^{S}\left(f\right)\ne\:{G}_{n,i}^{S}\left(f\right)\) or \(\:{G}_{m,i}^{B}\left(f\right)\ne\:{G}_{n,i}^{B}\left(f\right)\) ), violating the assumption underlying the standard TSR method. The violation means the site-specific response is not fully canceled and instead distorts the calculated path-effect term, which reduces the correlation between the surface-TSR and the borehole-TSR, directly influencing the interstation correlation of SBSR fluctuations (upper panel, magenta line). Within the same frequency bands (3–4 Hz and 6–9.5 Hz) where the path-effect correlation is low, the SBSR fluctuation correlation for MYGH04–MYGH13 becomes strongly negative, approaching − 0.5. This finding indicates an apparent inverse correlation between the SBSR fluctuations at the two sites. However, this inverse correlation should be interpreted with caution. It occurs precisely where the underlying assumption of the analysis, namely, that the common path effect is successfully isolated and removed, breaks down. Given that this negative correlation occurs precisely where the underlying assumption of the analysis—that the common path effect is successfully isolated—breaks down, we interpret this result as a methodological artifact rather than a true physical phenomenon. The inclusion of such strong negative correlations from a subset of station pairs explains the greater overall statistical dispersion observed in Dataset B, which is reflected by the wider whisker range in the upper panel of Fig. 5(b) compared with that of Fig. 5(a). Within the frequency bands where the path-effect correlation is high (i.e., below 3 Hz and between 4.5 and 6 Hz), the interstation SBSR fluctuation correlation is indeed close to zero. The results for Dataset A (Fig. 5(a)) provide a clearer picture than do those for Dataset B. For Dataset A, the path-effect correlation for MYGH04–MYGH13 (lower panel, magenta line) is consistently high, following the median trend. More importantly, the interstation correlation of SBSR fluctuations (upper panel, magenta line) remains consistently close to zero across the entire frequency range. Via the use of earthquake pairs where path-related differences are nearly eliminated, we isolate the site response to slight variations in the incoming wavefield. The fact that the resulting SBSR fluctuations at station MYGH13 are not correlated with those at the nearby MYGH04 station, even in frequency bands susceptible to notable azimuthal effects, demonstrates that this variability is a local site-specific effect. The analysis of the MYGH04–MYGH13 station pair provides two important observations. First, the negative correlation of SBSR fluctuations for Dataset B coincides with the frequency bands where the path-effect correlation is low. This suggests that the negative correlation is an artifact of applying the TSR method, in which the underlying assumptions are violated by notable, azimuth-dependent site effects. Second, the analysis of Dataset A, where path effects are minimal, reveals a near-zero correlation of SBSR fluctuations across all frequencies. These results indicate that the considerable azimuthal variability at MYGH13 is a localized characteristic attributable to the site conditions and is not spatially correlated with the nearby station. 5. Conclusion This study investigated the significant, azimuth-dependent variability of SBSRs using a novel TSR-based approach to distinguish between path- specific and site-specific effects. Our analysis of 293 KiK-net station pairs revealed a high correlation between path-effect terms derived from surface and borehole records but a near-zero interstation correlation of event-to-event SBSR fluctuations. The observed spatial decorrelation of SBSR variability provides strong evidence that its origin is predominantly local and site-specific. If this variability were driven by regional-scale path effects, a significant positive correlation would be expected between nearby stations. The absence of such a correlation indicates that the azimuthal dependence of site response is governed by site-specific complexities, such as 2D/3D geological structures. Furthermore, our experimental design, utilizing two distinct datasets, offers a potential framework for decomposing site-specific uncertainty. The variability observed using collocated earthquake pairs (Dataset A) can be interpreted as representing the aleatory uncertainty, the inherent randomness in the site response. In contrast, the variability from the different-azimuth pairs (Dataset B) encompasses this aleatory uncertainty plus the epistemic uncertainty associated with the systematic, unmodeled effects of azimuthal dependence. This suggests a pathway to validate the appropriateness and potential impact of applying non-ergodic site terms, by demonstrating how much of the total site-specific variability they can successfully explain. These results have direct implications for the next generation of seismic hazard assessment. By reclassifying a significant component of ground motion variability from random uncertainty to a deterministic, local site attribute, this work provides the physical basis and empirical data needed to develop more robust non-ergodic ground motion models (GMMs). Modern GMMs increasingly seek to break the ergodic assumption by separating the total variability into regional and site-specific components (Shahi and Baker, 2014 ; Baker and Bradley, 2014). This finding directly supports the ongoing development of non-ergodic GMMs. While these models already aim to isolate site-specific variability (Shahi and Baker, 2014 ), the indeterminate origin of its azimuthal component has been a key source of epistemic uncertainty. By demonstrating that this azimuthal dependence is a stable, site-specific phenomenon, our work provides a physical validation for treating this variability as a deterministic local attribute, paving the way for its explicit incorporation into future models. While this study provides a new framework for understanding site-specific variability, it is important to acknowledge its limitations. First, our methodology, while effective at identifying the origin of azimuthal dependence (site versus path), does not fully characterize the specific physical mechanisms (e.g., scattering from topography versus basin resonance) at each individual site. Second, our findings are based on the dense KiK-net in Japan, and their generalizability to other tectonic regions requires further investigation. Furthermore, our findings are based on the dense KiK-net in Japan; their generalizability to other tectonic regions requires further investigation. Future work should therefore focus on applying this methodology to other regions and integrating the results with detailed geological and topographical data to build more quantitative models of these 2D/3D effects. Such efforts are essential for developing the practical parameterizations needed to refine non-ergodic models (Lavrentiadis et al., 2023 ) and will ultimately lead to more reliable PSHA. Abbreviations 1D/2D/3D One-/two-/three-dimensional CSMIP California Strong-Motion Instrumentation Program DTR Dual twofold spectral ratio ETF Empirical transfer function \(\:F\) Fourier spectra of the observation records \(\:G\) Ground amplification function GMM Ground motion model IQR Interquartile range JMA Japan Meteorological Agency KiK-net Kiban–Kyoshin network ln Natural logarithm MANOVA Multivariate analysis of variance M J Magnitude defined by the JMA NIED National Research Institute for Earth Science and Disaster Resilience \(\:P\) Path transfer function PGA Peak ground acceleration PSHA Probabilistic seismic hazard analysis RSF Relative site factor \(\:S\) Source spectra SBSR Surface-to-borehole spectral ratio SD Standard deviation TSR Twofold spectral ratio TTF Theoretical transfer function Declarations Ethics approval and consent to participate Not applicable Consent for publication Not applicable Authors' contributions The author conceived the study, analyzed the data, and wrote the manuscript. Dr. Yamamoto, Dr. Tokumitsu, and Dr. Uchiyama contributed to the development of the research direction through regular discussions and provided critical feedback on the manuscript. Prof. Takai, as the research supervisor, provided overall guidance for this research and reviewed the manuscript. All the authors have read and approved the final manuscript. Acknowledgements We thank the NIED for providing the KiK-net observation records. Availability of data and materials The strong-motion data used in this study were obtained from the Kiban–Kyoshin network (KiK-net), operated by the National Research Institute for Earth Science and Disaster Resilience (NIED). The data are publicly available from the NIED website at https://www.kyoshin.bosai.go.jp/ . The analysis code used in this study is available from the corresponding author upon reasonable request. References Aki K, Richards PG (2002) Quantitative Seismology, 2nd edn. University Science Books Archuleta RJ, Seale SH, Sangas PV, Baker LM, Swain ST (1992) Garner Valley downhole array of accelerometers: Instrumentation and preliminary data analysis. Bull Seismol Soc Am 82:1592–1621. 10.1785/BSSA0820041592 Baker JW, Bradley BA (2017) Intensity Measure Correlations Observed in the NGA-West2 Database, and Dependence of Correlations on Rupture and Site Parameters. Earthq Spectra 33(1):145–156. 10.1193/060716eqs095m Bard PY, Bouchon M (1985) The two-dimensional resonance of sediment-filled valleys. Bull Seismol Soc Am 75(2):519–541. 10.1785/BSSA0750020519 Field EH, Jacob KH (1995) A comparison and test of various site-response estimation techniques, including three that are not reference-site dependent. Bull Seismol Soc Am 85:1127–1143. 10.1785/BSSA0850041127 Hair JF Jr, Black WC, Babin B, Anderson RE (2018) Multivariate Data Analysis, Eighth Edition. Cengage Learning EMEA Ikeura T (2020) A method for investigating the attenuation characteristics of strong motions without attenuation functions: Application to strong motions in the Tohoku and Kanto regions caused by M7 class plate boundary earthquakes. J Japan Association Earthq Eng. 10.5610/jaee.20.7_133 (in Japanese) 20:7_133-7_157 Izumitani Y (2000) Qs-value in southern Kyushu evaluated from double spectral ratio of strong motion records. Journal of Japan Society of Civil Engineers, 2000:225–230. 10.2208/jscej.2000.640_225 (in Japanese) Joyner WB, Warrick RE, Oliver AA III (1976) Analysis of seismograms from a downhole array in sediments near San Francisco Bay. Bull Seismol Soc Am 66(3):937–958. 10.1785/BSSA0660030937 Lavrentiadis G, Abrahamson NA, Kuehn NM (2023) A non-ergodic effective amplitude ground-motion model for California. Bull Earthq Eng 21:5233–5264. 10.1007/s10518-021-01206-w Matsuzawa T, Takeo M, Ide S, Iio Y, Ito H, Imanishi K, Horiuchi S (2004) Estimation of the S-wave attenuation in the western Nagano region from twofold spectral ratio. J Seismological Soc Japan Ser 2 56:75–88. 10.4294/zisin1948.56.1_75 (in Japanese) National Research Institute for Earth Science and Disaster Resilience (2019) NIED K-NET, KiK-net. 10.17598/NIED.0004 Noguchi S, Sasatani T (2011) Nonlinear soil response and its effects on strong ground motions during the 2003 Miyagi-Oki intraslab earthquake. J Seismological Soc Japan Ser 2 63:165–187. 10.4294/zisin.63.165 (in Japanese) Parker GA, Stewart JP, Boore DM, Atkinson GM, Hassani B (2021) NGA-subduction global ground motion models with regional adjustment factors. Earthq Spectra 38:456–493. 10.1177/87552930211034889 Pedregosa F, Varoquaux G, Gramfort A, Michel V, Thirion B, Grisel O, Blondel M, Prettenhofer P, Weiss R, Dubourg V, Vanderplas J, Passos A, Cournapeau D, Brucher M, Perrot M, Duchesnay É (2011) Scikit-learn: Machine Learning in Python. J Mach Learn Res 12:2825–2830 Pilz M, Cotton F (2019) Does the one-dimensional assumption hold for site response analysis? A study of seismic site responses and implication for ground motion assessment using KiK-Net strong-motion data. Earthq Spectra 35:883–905. 10.1193/050718EQS113M Seabold S, Perktold J (2010) statsmodels: Econometric and statistical modeling with python. In: Proceedings of the 9th Python in Science Conference Sgobba S, Felicetta C, Bortolotti T, Menafoglio A, Lanzano G, Pacor F (2024) A geostatistical modelling of empirical amplification functions and related site proxies for shaking scenarios in central Italy. Soil Dyn Earthq Eng 179:108496. 10.1016/j.soildyn.2024.108496 Shahi SK, Baker JW (2014) NGA-West2 Models for Ground Motion Directionality. Earthq Spectra 30(3):1285–1300. 10.1193/040913EQS097M Shiomi K (2012) An improved method for estimating the installation azimuth information of NIED Hi-net borehole seismometers. Rep Natl Res Inst Earth Sci Disaster Resil 80:1–20 (in Japanese) Shima E (1962) Modification of Seismic Waves in Superficial Layers as Verified by Comparative Observations on and Beneath the Surface. Bull Earthq Res Inst Tokyo 40:187–259 Tao Y, Rathje E (2020) The importance of distinguishing pseudoresonances and outcrop resonances in downhole array data. Bull Seismol Soc Am 110:288–294. 10.1785/0120190097 Thompson EM, Baise LG, Tanaka Y, Kayen RE (2012) A taxonomy of site response complexity. Soil Dyn Earthq Eng 41:32–43. 10.1016/j.soildyn.2012.04.005 Ueno H, Hatakeyama S, Aketagawa T, Funazaki A, Hamada N (2002) Improvement of hypocenter determination procedures in the Japan Meteorological Agency. Seismological Bull Japan Meteorological Agency 65:123–134 (in Japanese) Supplementary Files 0GraphicalAbstract.png Cite Share Download PDF Status: Posted Version 1 posted You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-7806012","acceptedTermsAndConditions":true,"allowDirectSubmit":true,"archivedVersions":[],"articleType":"Research Article","associatedPublications":[],"authors":[{"id":536676581,"identity":"8c05f830-a1e9-49dd-9d2b-55c86f9aa601","order_by":0,"name":"Masashi Aoki","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAABCElEQVRIiWNgGAWjYFACNhAhwcMPZfEwMDA2gBjMBLTYyEk2gFgJxGtJMzY4ANFCGPC3t6Vu/PHncOLm2z1mDxh/2MnItzc3MPyoYGA3x6FF4syxY7d52w4nbrtzxtyAISGZx+DMwQbGnjMMzJYNOPTcSG+7zdgA1HIjx0yCIYGZx0AisYGZsY2BGehUrEAeqOUm2GEzwFrqeeTnP8SvxeBG2rEbPGxA70uAtRzmYbjBiF+L4ZljaUC/2MhJ3DlWbpCQdhzol8SGgz1nJHD6Re54mxnQYcConN287cEHm2p7+fbjDx/8qLBJxhViCCDBgIgWoJMkkg2I0oIM7AhrGQWjYBSMghECAJIdW06zFauNAAAAAElFTkSuQmCC","orcid":"","institution":"Taisei Corporation","correspondingAuthor":true,"prefix":"","firstName":"Masashi","middleName":"","lastName":"Aoki","suffix":""},{"id":536676582,"identity":"0a3bba23-7c0d-4d9f-b1ad-da7dcdd8fbe5","order_by":1,"name":"Yu Yamamoto","email":"","orcid":"","institution":"Taisei Corporation","correspondingAuthor":false,"prefix":"","firstName":"Yu","middleName":"","lastName":"Yamamoto","suffix":""},{"id":536676583,"identity":"ca07f99f-b7c9-4619-a0f3-52930aaf0d46","order_by":2,"name":"Ryoichi Tokumitsu","email":"","orcid":"","institution":"Taisei Corporation","correspondingAuthor":false,"prefix":"","firstName":"Ryoichi","middleName":"","lastName":"Tokumitsu","suffix":""},{"id":536676584,"identity":"bf5e562c-a175-413f-ad08-470165b716fa","order_by":3,"name":"Yasuo Uchiyama","email":"","orcid":"","institution":"Taisei Corporation","correspondingAuthor":false,"prefix":"","firstName":"Yasuo","middleName":"","lastName":"Uchiyama","suffix":""},{"id":536676585,"identity":"2a699193-ab99-44e4-8490-79eab8f67068","order_by":4,"name":"Nobuo Takai","email":"","orcid":"","institution":"Hokkaido University: Hokkaido Daigaku","correspondingAuthor":false,"prefix":"","firstName":"Nobuo","middleName":"","lastName":"Takai","suffix":""}],"badges":[],"createdAt":"2025-10-08 09:03:37","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-7806012/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-7806012/v1","draftVersion":[],"editorialEvents":[],"editorialNote":"","failedWorkflow":false,"files":[{"id":95355573,"identity":"6bfdecc9-3f73-4a99-b06b-a979928cd6ed","added_by":"auto","created_at":"2025-11-07 06:24:03","extension":"xml","order_by":6,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":10413,"visible":true,"origin":"","legend":"","description":"","filename":"epspEPSPD2500319.xml","url":"https://assets-eu.researchsquare.com/files/rs-7806012/v1/1fc2cea26622984c3fd8aaeb.xml"},{"id":95355563,"identity":"b7ca53f3-9f16-4aab-b0f6-92bf1e94a471","added_by":"auto","created_at":"2025-11-07 06:24:02","extension":"xml","order_by":7,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":1195,"visible":true,"origin":"","legend":"","description":"","filename":"EPSPD250031911401.go.xml","url":"https://assets-eu.researchsquare.com/files/rs-7806012/v1/8487d54344f2fd8484ad5947.xml"},{"id":95355566,"identity":"bf75fff3-fc95-4d5c-9c6c-11459330a43b","added_by":"auto","created_at":"2025-11-07 06:24:03","extension":"xml","order_by":8,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":917,"visible":true,"origin":"","legend":"","description":"","filename":"EPSPD2500319Import.xml","url":"https://assets-eu.researchsquare.com/files/rs-7806012/v1/fb2862e52ee3c45de45cbea9.xml"},{"id":95524398,"identity":"2adae031-2c64-4e91-b02a-d85940a83d82","added_by":"auto","created_at":"2025-11-10 10:02:43","extension":"xml","order_by":10,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":90043,"visible":true,"origin":"","legend":"","description":"","filename":"EPSPD25003191enriched.xml","url":"https://assets-eu.researchsquare.com/files/rs-7806012/v1/66c37a7babd2e4e4052ed699.xml"},{"id":95355567,"identity":"33c2d224-f784-492a-9747-5789b1b632aa","added_by":"auto","created_at":"2025-11-07 06:24:03","extension":"png","order_by":11,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":354865,"visible":true,"origin":"","legend":"","description":"","filename":"1histograms.png","url":"https://assets-eu.researchsquare.com/files/rs-7806012/v1/4952110e7544b100b8fdc320.png"},{"id":95355568,"identity":"ee360833-7db1-4941-93bc-40b4e5e35f0f","added_by":"auto","created_at":"2025-11-07 06:24:03","extension":"png","order_by":12,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":1599790,"visible":true,"origin":"","legend":"","description":"","filename":"2PreliminaryAnalysiis.png","url":"https://assets-eu.researchsquare.com/files/rs-7806012/v1/286407ed7f9c7c298ace260f.png"},{"id":95355575,"identity":"72134a78-dfa3-4b19-b9af-b462fec8595b","added_by":"auto","created_at":"2025-11-07 06:24:03","extension":"png","order_by":13,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":188903,"visible":true,"origin":"","legend":"","description":"","filename":"3schematicdiagram.png","url":"https://assets-eu.researchsquare.com/files/rs-7806012/v1/ffef38400a58c451cb3bac0e.png"},{"id":95355564,"identity":"7b5289aa-9ced-461b-98ca-934db7e2431b","added_by":"auto","created_at":"2025-11-07 06:24:02","extension":"png","order_by":14,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":585557,"visible":true,"origin":"","legend":"","description":"","filename":"4spacialdistribution.png","url":"https://assets-eu.researchsquare.com/files/rs-7806012/v1/0d22050b1476c7fb0dbb602d.png"},{"id":95355582,"identity":"c1325a3f-c7e1-4a1b-9a8d-e20b6dcbe71c","added_by":"auto","created_at":"2025-11-07 06:24:04","extension":"png","order_by":15,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":713938,"visible":true,"origin":"","legend":"","description":"","filename":"5correlation.png","url":"https://assets-eu.researchsquare.com/files/rs-7806012/v1/1cc52fa4d23d5979fb2acfc8.png"},{"id":95355569,"identity":"fc853d71-1507-48f8-9997-213e9f7f238f","added_by":"auto","created_at":"2025-11-07 06:24:03","extension":"png","order_by":16,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":73708,"visible":true,"origin":"","legend":"","description":"","filename":"Online1histograms.png","url":"https://assets-eu.researchsquare.com/files/rs-7806012/v1/c3a81df2e1db599a0ba84460.png"},{"id":95525251,"identity":"119fec85-02ad-4f8b-aa56-362fef155ea7","added_by":"auto","created_at":"2025-11-10 10:04:37","extension":"png","order_by":17,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":229432,"visible":true,"origin":"","legend":"","description":"","filename":"Online2PreliminaryAnalysiis.png","url":"https://assets-eu.researchsquare.com/files/rs-7806012/v1/3e8444a82eb4360944c1e6ec.png"},{"id":95355570,"identity":"7e0dc005-23b4-4a54-8d6f-96d5a911763a","added_by":"auto","created_at":"2025-11-07 06:24:03","extension":"png","order_by":18,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":34705,"visible":true,"origin":"","legend":"","description":"","filename":"Online3schematicdiagram.png","url":"https://assets-eu.researchsquare.com/files/rs-7806012/v1/7a625672f8c8ec389c3a2317.png"},{"id":95355577,"identity":"36523dc4-1a3b-4cbc-8d2b-d0b4e93ac189","added_by":"auto","created_at":"2025-11-07 06:24:03","extension":"png","order_by":19,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":104278,"visible":true,"origin":"","legend":"","description":"","filename":"Online4spacialdistribution.png","url":"https://assets-eu.researchsquare.com/files/rs-7806012/v1/e5453c40b620e8578c8aa05f.png"},{"id":95355562,"identity":"c71c283f-0e74-4f48-b641-b7424713a185","added_by":"auto","created_at":"2025-11-07 06:24:02","extension":"png","order_by":20,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":77514,"visible":true,"origin":"","legend":"","description":"","filename":"Online5correlation.png","url":"https://assets-eu.researchsquare.com/files/rs-7806012/v1/b5d98d46d771ec986bec085f.png"},{"id":95355578,"identity":"7724d17d-38e2-49e3-be63-bf80545321a3","added_by":"auto","created_at":"2025-11-07 06:24:03","extension":"xml","order_by":21,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":88609,"visible":true,"origin":"","legend":"","description":"","filename":"EPSPD25003191structuring.xml","url":"https://assets-eu.researchsquare.com/files/rs-7806012/v1/2ada777493d5ea347173784d.xml"},{"id":95355579,"identity":"b240c93a-67b0-45b9-84a3-44a9650bbd53","added_by":"auto","created_at":"2025-11-07 06:24:04","extension":"html","order_by":22,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":98828,"visible":true,"origin":"","legend":"","description":"","filename":"earlyproof.html","url":"https://assets-eu.researchsquare.com/files/rs-7806012/v1/e996d3486be8145a4ce03bb0.html"},{"id":95355576,"identity":"d29a11f9-6b42-4a36-87cd-c77ad5b6b5c8","added_by":"auto","created_at":"2025-11-07 06:24:03","extension":"png","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":354865,"visible":true,"origin":"","legend":"\u003cp\u003eDistributions of key parameters for the selected dataset. The panels show histograms of (a) JMA-defined magnitude (M\u003csub\u003eJ\u003c/sub\u003e), (b) focal depth, (c) surface peak ground acceleration (PGA), (d) epicentral distance, (e) hypocentral distance, and (f) geometric incidence angle for the 152,352 records used in this study.\u003c/p\u003e","description":"","filename":"1histograms.png","url":"https://assets-eu.researchsquare.com/files/rs-7806012/v1/ef159062718a0f9bf05c77e4.png"},{"id":95355561,"identity":"37f99ba9-1dcf-4918-8b3f-05a1dc2eb659","added_by":"auto","created_at":"2025-11-07 06:24:02","extension":"png","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":1599790,"visible":true,"origin":"","legend":"\u003cp\u003ePredominance of azimuthal effects on surface-to-borehole spectral ratio (SBSR) variability. (a) Frequency-dependent standard deviation of ln(SBSR) for all 693 stations. (b) Mean SBSRs for station MYGH13, grouped by back-azimuth. (c) Example of SBSR peak clustering using a Gaussian mixture model. (d) Scatter plot of Pillai’s trace statistics from a MANOVA, comparing the effects of back-azimuth and geometric incidence angle on the lowest-frequency SBSR peak cluster for 347 stations. The MANOVA quantitatively confirms that back-azimuth has a more notable influence.\u003c/p\u003e","description":"","filename":"2PreliminaryAnalysiis.png","url":"https://assets-eu.researchsquare.com/files/rs-7806012/v1/6aa432f3be829358253d5ac8.png"},{"id":95355581,"identity":"9a2d5457-d185-43a2-885c-9a2f85516811","added_by":"auto","created_at":"2025-11-07 06:24:04","extension":"png","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":188903,"visible":true,"origin":"","legend":"\u003cp\u003eSchematic illustration of the modified twofold spectral ratio (TSR) framework for isolating site-specific variability. The ratio of the surface-TSR to the borehole-TSR is designed to cancel the common differential path-effect term. This leaves the ratio of the inter-event surface-to-borehole spectral ratio (SBSR) fluctuations between stations i and j, allowing for a direct evaluation of the spatial coherence of site-specific variability.\u003c/p\u003e","description":"","filename":"3schematicdiagram.png","url":"https://assets-eu.researchsquare.com/files/rs-7806012/v1/5a08cc2ddb4a4619ea49b75a.png"},{"id":95355580,"identity":"540e1d73-2edc-4b11-b319-60a1d85fc188","added_by":"auto","created_at":"2025-11-07 06:24:04","extension":"png","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":585557,"visible":true,"origin":"","legend":"\u003cp\u003eSpatial distribution of the selected station pairs and earthquake datasets. (a) Map showing the 293 station pairs selected for analysis, with an interstationdistance smaller than 20 km. (b) Example distribution of collocatedearthquake pairs (Dataset A) for the MYGH04–MYGH13 station pair. (c) Example distribution of common-distance, different-path earthquake pairs (Dataset B) for the same station pair.\u003c/p\u003e","description":"","filename":"4spacialdistribution.png","url":"https://assets-eu.researchsquare.com/files/rs-7806012/v1/d943f0210e47287be0f2ef5e.png"},{"id":95525954,"identity":"10631f25-4839-477c-a34a-377219ab01d1","added_by":"auto","created_at":"2025-11-10 10:05:57","extension":"png","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":713938,"visible":true,"origin":"","legend":"\u003cp\u003eNear-zero interstation correlation of surface-to-borehole spectral ratio (SBSR) fluctuations indicates a site-specific origin for azimuthal variability. The panels show frequency-dependent Pearson’s correlation coefficients for (a) collocated earthquake pairs (Dataset A) and (b) common-distance, different-path pairs (Dataset B). The upper subpanels show the interstation correlation of SBSR fluctuations. The lower subpanels show the path-effect correlation between the surface-twofold spectral ratio (TSR) and borehole-TSR. The solid line indicates the median, the dark shaded area denotes the interquartile range (IQR), and the light shaded area represents 1.5×IQR. The results for the MYGH04–MYGH13 station pair are highlighted by magenta lines.\u003c/p\u003e","description":"","filename":"5correlation.png","url":"https://assets-eu.researchsquare.com/files/rs-7806012/v1/6f9b4ffeada9de75719ea1c5.png"},{"id":97667642,"identity":"c5d60e49-923e-441a-af17-af9a680b2e29","added_by":"auto","created_at":"2025-12-08 09:23:57","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":3628633,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-7806012/v1/91a12bce-bb2d-40bc-9222-146d95c081ae.pdf"},{"id":95355571,"identity":"e9ec7d18-d6a4-4122-8561-ae635f81ae7f","added_by":"auto","created_at":"2025-11-07 06:24:03","extension":"png","order_by":9,"title":"","display":"","copyAsset":false,"role":"supplement","size":150338,"visible":true,"origin":"","legend":"","description":"","filename":"0GraphicalAbstract.png","url":"https://assets-eu.researchsquare.com/files/rs-7806012/v1/271eae81366444ff0a2e8115.png"}],"financialInterests":"","formattedTitle":"Investigating the Origin of Azimuthal Site Response Variability using KiK-net Data: A Twofold Spectral Ratio Approach","fulltext":[{"header":"1. Introduction","content":"\u003cp\u003eSite-specific ground motion amplification is a critical factor in earthquake-resistant design and seismic hazard assessment. One-dimensional (1D) site response analysis has long been a fundamental tool in geotechnical earthquake engineering practice, simplifying wave propagation to vertically propagating shear waves through horizontal soil layers (Kramer, 1996). Within this framework, site amplification is characterized by a theoretical transfer function (TTF). The TTF assumes ergodicity, where a site\u0026rsquo;s response is expected to be stable and time-invariant, representing a unique and predictable property of its local geology.\u003c/p\u003e\u003cp\u003eA key metric for empirically measuring the TTF is the surface-to-borehole spectral ratio (SBSR). Pioneering studies demonstrated that by analyzing simultaneous recordings from surface and subsurface instruments, it was possible to compute actual amplification functions\u0026mdash;often treated as the empirical transfer function (ETF) of the site\u0026mdash;and validate the core principles of site response (Shima, \u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e1962\u003c/span\u003e; Joyner et al., \u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e1976\u003c/span\u003e). Extensive vertical array networks, such as the Kiban\u0026ndash;Kyoshin network (KiK-net) in Japan, the California Strong-Motion Instrumentation Program (CSMIP), and the Strong-Motion Downhole Array in the Taipei Basin, have provided a substantial volume of data, facilitating the widespread application of the SBSR method.\u003c/p\u003e\u003cp\u003eNevertheless, the interpretation of the SBSR as a pure transfer function is complicated. Ideally, the SBSR would capture a true resonance of the site, which corresponds to the amplification at its natural frequencies under outcrop conditions. However, even within an idealized 1D framework, the interference from downward-propagating waves creates spectral troughs in the borehole recording, which can result in pseudoresonance peaks in the calculated SBSR. In practice, this theoretical complexity is further compounded, as observed borehole seismograms are also contaminated by complex wavefields, such as scattered waves from 3D geological irregularities (Thompson et al., \u003cspan citationid=\"CR23\" class=\"CitationRef\"\u003e2012\u003c/span\u003e; Tao and Rathje, \u003cspan citationid=\"CR22\" class=\"CitationRef\"\u003e2020\u003c/span\u003e). These physical complexities result in significant event-to-event variability, a phenomenon that directly challenges the ergodic assumption of 1D models, even for weak motions within the linear response range (Archuleta et al., \u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e1992\u003c/span\u003e; Pilz and Cotton, \u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e2019\u003c/span\u003e). The overall variability reflects the limitations of 1D models, contributing to epistemic uncertainty in ground motion models (GMMs) and probabilistic seismic hazard analysis (PSHA) (Parker et al., \u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e2021\u003c/span\u003e).\u003c/p\u003e\u003cp\u003eThe event-to-event variability in site responses can be attributed to physical factors including the wave type, back-azimuth, and incidence angle (Field and Jacob, \u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e1995\u003c/span\u003e). The presence of azimuth-dependent variability is a critical issue for modern data-driven methods that aim to break the ergodic assumption by separating local site amplification from regional path effects. This gives rise to a fundamental ambiguity, as the physical origin of this azimuthal dependence\u0026mdash;whether it is attributable to the site term or the path term\u0026mdash;has remained unresolved. Quantitatively distinguishing between these two sources is a crucial step for advancing site characterization. The practical importance of this distinction is evident in the development of non-ergodic GMMs. For instance, Ikeura (\u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e2020\u003c/span\u003e) proposed a method to investigate path characteristics without assuming a uniform attenuation function by first removing site effects, which were quantified as relative site factors (RSFs). The stability and validity of such an RSF, however, is directly challenged if it exhibits an azimuthal dependence of unknown origin. This highlights the need to resolve whether the dependence is a site-related or path-related phenomenon.\u003c/p\u003e\u003cp\u003eThe objective of this study is to investigate the physical origin of the observed azimuthal dependence of SBSRs. We employ an analytical approach based on the twofold spectral ratio (TSR) method, adapted for pairs of nearby vertical array stations, to systematically decompose site response variability into path-related and site-related components.\u003c/p\u003e"},{"header":"2. Data and Preliminary Analysis","content":"\u003cp\u003eThe dataset and preprocessing steps are described in this section. It also presents a preliminary analysis to characterize SBSR variability and to identify the primary factors influencing this variability.\u003c/p\u003e\u003cdiv id=\"Sec3\" class=\"Section2\"\u003e\u003ch2\u003e2.1. KiK-net Strong-Motion Data\u003c/h2\u003e\u003cp\u003eData were obtained from the KiK-net, a nationwide strong-motion seismograph network in Japan operated by the National Research Institute for Earth Science and Disaster Resilience (NIED). The network comprises approximately 700 stations, each with surface and borehole three-component accelerometers. This configuration is ideal for directly observing site amplification effects and provides an extensive dataset for our research (NIED 2019).\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec4\" class=\"Section2\"\u003e\u003ch2\u003e2.2. Earthquake Record Selection\u003c/h2\u003e\u003cp\u003eTo compile a comprehensive, high-quality dataset, we selected earthquake records on the basis of several criteria, with the primary objective of analyzing the linear or quasilinear behavior of site responses. The selection criteria are as follows:\u003c/p\u003e\u003cp\u003e\u003col\u003e\u003cspan\u003e\u003cli\u003e\u003cp\u003e\u003cb\u003eEarthquake Magnitude\u003c/b\u003e: To ensure a sufficient signal-to-noise ratio for stable spectral analysis, only earthquakes with a Japan Meteorological Agency (JMA) -defined magnitude (M\u003csub\u003eJ\u003c/sub\u003e) of 4.0 or greater were included.\u003c/p\u003e\u003c/li\u003e\u003c/span\u003e\u003cspan\u003e\u003cli\u003e\u003cp\u003e\u003cb\u003ePeak Ground Acceleration (PGA)\u003c/b\u003e: Records were selected if the surface PGA varied between 1.0 and 50.0 cm/s\u0026sup2;. The upper limit was intended to exclude records that may have been significantly affected by soil nonlinearity, following the approach of Noguchi and Sasatani (\u003cspan citationid=\"CR13\" class=\"CitationRef\"\u003e2011\u003c/span\u003e).\u003c/p\u003e\u003c/li\u003e\u003c/span\u003e\u003cspan\u003e\u003cli\u003e\u003cp\u003e\u003cb\u003eEvent Uniqueness\u003c/b\u003e: To isolate the ground motion from a single source, records were screened to exclude instances in which multiple earthquakes occurred in close succession by examining theoretical travel-time predictions via the JMA2001 travel-time table (Ueno et al. \u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e2002\u003c/span\u003e).\u003c/p\u003e\u003c/li\u003e\u003c/span\u003e\u003c/ol\u003e\u003c/p\u003e\u003cp\u003eAfter applying these criteria, our final dataset comprised 152,352 records from 6,973 earthquakes observed at 693 stations. The distribution of key dataset parameters is shown in Fig.\u0026nbsp;1, which includes histograms of M\u003csub\u003eJ\u003c/sub\u003e, focal depth, surface PGA, epicentral distance, hypocentral distance, and geometric incidence angle. The geometric incidence angle is the angle calculated from the straight-line path between the hypocenter and the station, ignoring reflections and refractions at layer boundaries.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec5\" class=\"Section2\"\u003e\u003ch2\u003e2.3. Data Preprocessing\u003c/h2\u003e\u003cp\u003eAll acceleration time histories were subjected to a series of standard preprocessing steps. Baseline correction was performed by subtracting the mean of the entire record. A critical step in analyzing borehole records is correcting the orientation of horizontal components. In this study, we applied correction angle data published by Shiomi (\u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e2012\u003c/span\u003e). For the 33 stations lacking such data, an independent orientation correction was performed by maximizing the coherence between the surface and borehole horizontal recordings in a low-frequency band.\u003c/p\u003e\u003cp\u003eTo analyze the SH-wave response, the two horizontal components were rotated to the great-circle path from the epicenter to the station to obtain the transverse component. The SBSR was then computed by dividing the Fourier amplitude spectrum of the transverse-component record at the surface by that in the borehole. To obtain a stable spectral ratio and clarify the trend of site response, the resulting ratio was smoothed using a Parzen window (Noguchi and Sasatani, \u003cspan citationid=\"CR13\" class=\"CitationRef\"\u003e2011\u003c/span\u003e). A bandwidth of 1.0 Hz was selected to enhance the stability of the spectral ratio across the frequency range of interest.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec6\" class=\"Section2\"\u003e\u003ch2\u003e2.4. Preliminary Analysis of SBSR Variability\u003c/h2\u003e\u003cp\u003eTo understand the nature of the observed site response, the preliminary analysis focused on the 1\u0026ndash;10 Hz frequency range, which is a band of primary interest for many engineering applications (Kramer, 1996). The analysis results are shown in Fig.\u0026nbsp;2. The standard deviation (SD) of the natural logarithm of the SBSR (ln(SBSR)) for all 693 stations is shown in Fig.\u0026nbsp;2(a). The results for the station with the largest maximum SD (MYGH13) are highlighted in red. This panel reveals that the SBSR variability is generally small at low frequencies but increases significantly with frequency. The azimuthal dependence of the response is shown in Fig.\u0026nbsp;2(b), which shows the mean SBSR for station MYGH13, grouped into four back-azimuth quadrants. The spectral shape clearly differs for each direction.\u003c/p\u003e\u003cp\u003eTo objectively classify the observed SBSR peaks, a Gaussian mixture model\u0026mdash;a probabilistic model used to represent subpopulations within an overall population\u0026mdash;was applied using the scikit-learn library in Python (Pedregosa et al, \u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e2011\u003c/span\u003e) (Fig.\u0026nbsp;2(c)). The number of clusters in the model was defined on the basis of the average number of peaks observed per event at a given station. In the subsequent statistical analysis, we focused on the cluster with the lowest mean frequency. The variability in the lowest-frequency cluster, which could be caused by either a true resonance or pseudoresonance peak resulting from downward-propagating waves (Thompson et al. \u003cspan citationid=\"CR23\" class=\"CitationRef\"\u003e2012\u003c/span\u003e; Tao and Rathje \u003cspan citationid=\"CR22\" class=\"CitationRef\"\u003e2020\u003c/span\u003e), was considered the fundamental feature in this analysis. The clustering results for a representative station are shown in Fig.\u0026nbsp;2(c).\u003c/p\u003e\u003cp\u003eFinally, one-way multivariate analysis of variance (MANOVA) of the peaks belonging to the lowest-frequency cluster was conducted to quantitatively compare the influences of the back-azimuth and geometric incidence angle, implemented with the statsmodels library in Python (Seabold and Perktold, \u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e2010\u003c/span\u003e). The results are shown in the scatter plot of Pillai\u0026rsquo;s trace statistics in Fig.\u0026nbsp;2(d). MANOVA is a statistical method for comparing the means of two or more populations and for determining whether these differences are statistically significant. This technique has also been applied in fields such as geostatistical modeling (Sgobba et al. \u003cspan citationid=\"CR18\" class=\"CitationRef\"\u003e2024\u003c/span\u003e). Pillai's trace was selected as the test statistic for our MANOVA because it is particularly robust against violations of assumptions, such as unequal sample sizes between groups, making it suitable for analyzing earthquake datasets (Hair et al., \u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e2018\u003c/span\u003e). Among the 693 total stations, 347 stations provided sufficient data for a valid analysis. The results revealed that the statistics for the back-azimuth are more notable than those for the geometric incidence angle at most stations. Specifically, the back-azimuth was the dominant factor at 303 stations (87.3%), whereas the geometric incidence angle dominated at only 44 (12.7%).\u003c/p\u003e\u003cp\u003eThis preliminary analysis empirically establishes that the back-azimuth is the primary factor controlling the observed SBSR variability, raising the critical question of whether this considerable azimuthal dependence is caused by complexities along the wave propagation path or by local site conditions.\u003c/p\u003e\u003c/div\u003e"},{"header":"3. Methodology: Modified Twofold Spectral Ratio Approach","content":"\u003cp\u003eThe preliminary analysis in the previous section established that the observed SBSR variability is controlled primarily by the back-azimuth of incoming seismic waves, necessitating a distinction between two potential sources for the observed azimuthal dependence, namely, path- and site-specific effects. The methodology presented herein is therefore designed to quantitatively separate these two effects for analysis.\u003c/p\u003e\u003cdiv id=\"Sec8\" class=\"Section2\"\u003e\u003ch2\u003e3.1. Analytical Framework\u003c/h2\u003e\u003cp\u003eThe TSR method was established to estimate path-dependent seismic wave characteristics, most notably the quality factor (Q value) for anelastic attenuation, by systematically canceling common source and site amplification effects from observed ground motion spectra (Matsuzawa et al. \u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e2004\u003c/span\u003e; Izutani 2000). In this study, we adapted the TSR framework not to estimate path effects but to first isolate and then eliminate them, enabling a focused evaluation of site-specific contributions to SBSR variability.\u003c/p\u003e\u003cp\u003eThe Fourier spectrum of an observed ground motion record is typically modeled in the frequency domain as a convolution of source, path, and site effects. For a simplified point-source and far-field approximation, this can be expressed as a product. The observed Fourier spectra at the surface (\u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{F}_{m,i}^{S}\\)\u003c/span\u003e\u003c/span\u003e) and in the borehole (\u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{F}_{m,i}^{B}\\)\u003c/span\u003e\u003c/span\u003e) for station \u003cem\u003ei\u003c/em\u003e and earthquake event \u003cem\u003em\u003c/em\u003e are thus defined as follows (Aki and Richards, \u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e2002\u003c/span\u003e):\u003cdiv id=\"Equ1\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equ1\" name=\"EquationSource\"\u003e\n$$\\:{F}_{m,i}^{S}\\left(f\\right)={S}_{m}\\left(f\\right)\\bullet\\:{P}_{m,i}\\left(f\\right)\\bullet\\:{G}_{m,i}^{S}\\left(f\\right)$$\u003c/div\u003e\u003cdiv class=\"EquationNumber\"\u003e1\u003c/div\u003e\u003c/div\u003e\u003cdiv id=\"Equ2\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equ2\" name=\"EquationSource\"\u003e\n$$\\:{F}_{m,i}^{B}\\left(f\\right)={S}_{m}\\left(f\\right)\\bullet\\:{P}_{m,i}\\left(f\\right)\\bullet\\:{G}_{m,i}^{B}\\left(f\\right)$$\u003c/div\u003e\u003cdiv class=\"EquationNumber\"\u003e2\u003c/div\u003e\u003c/div\u003e\u003c/p\u003e\u003cp\u003ewhere \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{S}_{m}\\left(f\\right)\\)\u003c/span\u003e\u003c/span\u003e denotes the source spectrum for event \u003cem\u003em\u003c/em\u003e, \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{P}_{m,i}\\left(f\\right)\\)\u003c/span\u003e\u003c/span\u003e denotes the path transfer function from source \u003cem\u003em\u003c/em\u003e to station site \u003cem\u003ei\u003c/em\u003e, \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{G}_{m,i}^{S}\\left(f\\right)\\)\u003c/span\u003e\u003c/span\u003e is the surface site amplification factor at station \u003cem\u003ei\u003c/em\u003e, and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{G}_{m,i}^{B}\\left(f\\right)\\)\u003c/span\u003e\u003c/span\u003e is the corresponding site response term at the borehole depth. It is essential to recognize that the borehole recording does not constitute pure input motion. The wavefield at the borehole sensor is invariably contaminated by downward-propagating waves reflected from the free surface and any overlying low-velocity layers. This downgoing wave effect complicates the interpretation of borehole records for use as a reference for incident motion. The interference between upward- and downward-propagating waves can produce spectral troughs in borehole recordings, which in turn are manifested as apparent or pseudoresonance peaks in the calculated SBSR that are not representative of simple 1D resonance of the soil column (Thompson et al. \u003cspan citationid=\"CR23\" class=\"CitationRef\"\u003e2012\u003c/span\u003e; Tao and Rathje \u003cspan citationid=\"CR22\" class=\"CitationRef\"\u003e2020\u003c/span\u003e). The SBSR, which is the direct ETF between the surface and borehole, can be calculated as follows:\u003cdiv id=\"Equ3\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equ3\" name=\"EquationSource\"\u003e\n$$\\:{SBSR}_{m,i}\\left(f\\right)={F}_{m,i}^{S}\\left(f\\right)/{F}_{m,i}^{B}\\left(f\\right)={G}_{m,i}^{S}\\left(f\\right)/{G}_{m,i}^{B}\\left(f\\right)$$\u003c/div\u003e\u003cdiv class=\"EquationNumber\"\u003e3\u003c/div\u003e\u003c/div\u003e\u003c/p\u003e\u003cp\u003eIn this ratio, the shared source and path terms are theoretically eliminated, leaving a relatively stable characteristic of the vertical site structure, with a theoretical dependence on the incidence angle that is negligible for near-vertically propagating waves.\u003c/p\u003e\u003cp\u003eThe TSR method is applied to a pair of stations (\u003cem\u003ei\u003c/em\u003e and \u003cem\u003ej\u003c/em\u003e) and a pair of events (\u003cem\u003em\u003c/em\u003e and \u003cem\u003en\u003c/em\u003e). The TSR for surface recordings (\u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{TSR}^{S}\\)\u003c/span\u003e\u003c/span\u003e) is formulated by computing a ratio of spectral ratios. This operation cancels the source terms and, under the assumption of linear and event-invariant site responses, cancels the site terms, isolating the differential path effect:\u003cdiv id=\"Equ4\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equ4\" name=\"EquationSource\"\u003e\n$$\\:{TSR}^{S}\\left(f\\right)=\\frac{{F}_{m,i}^{S}\\left(f\\right)/{F}_{n,i}^{S}\\left(f\\right)}{{F}_{m,j}^{S}\\left(f\\right)/{F}_{n,j}^{S}\\left(f\\right)}=\\frac{{P}_{m,i}\\left(f\\right)/{P}_{n,i}\\left(f\\right)}{{P}_{m,j}\\left(f\\right)/{P}_{n,j}\\left(f\\right)}$$\u003c/div\u003e\u003cdiv class=\"EquationNumber\"\u003e4\u003c/div\u003e\u003c/div\u003e\u003c/p\u003e\u003cp\u003eThe same operation is applied to borehole records (\u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{TSR}^{B}\\)\u003c/span\u003e\u003c/span\u003e) to isolate the differential path-effect term:\u003cdiv id=\"Equ5\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equ5\" name=\"EquationSource\"\u003e\n$$\\:{TSR}^{B}\\left(f\\right)=\\frac{{F}_{m,i}^{B}\\left(f\\right)/{F}_{n,i}^{B}\\left(f\\right)}{{F}_{m,j}^{B}\\left(f\\right)/{F}_{n,j}^{B}\\left(f\\right)}=\\frac{{P}_{m,i}\\left(f\\right)/{P}_{n,i}\\left(f\\right)}{{P}_{m,j}\\left(f\\right)/{P}_{n,j}\\left(f\\right)}$$\u003c/div\u003e\u003cdiv class=\"EquationNumber\"\u003e5\u003c/div\u003e\u003c/div\u003e\u003c/p\u003e\u003cp\u003eThe central proposition of this study is that by computing the ratio of these two TSRs, i.e., the dual-TSR (DTR), common differential path effects, which both ratios are designed to capture, can be canceled.\u003cdiv id=\"Equ6\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equ6\" name=\"EquationSource\"\u003e\n$$\\:DTR\\left(f\\right)={TSR}^{S}\\left(f\\right)/{TSR}^{B}\\left(f\\right)$$\u003c/div\u003e\u003cdiv class=\"EquationNumber\"\u003e6\u003c/div\u003e\u003c/div\u003e\u003c/p\u003e\u003cp\u003eAs a direct result of this formulation, the resulting ratio is mathematically sensitive to any violation of the assumption that the empirical site response (i.e., the SBSR) is independent of the earthquake event. By substituting the definition of the SBSR in Eq.\u0026nbsp;(\u003cspan refid=\"Equ3\" class=\"InternalRef\"\u003e3\u003c/span\u003e) into Eq.\u0026nbsp;(\u003cspan refid=\"Equ6\" class=\"InternalRef\"\u003e6\u003c/span\u003e), the DTR can be rearranged into the following equivalent form:\u003cdiv id=\"Equ7\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equ7\" name=\"EquationSource\"\u003e\n$$\\:DTR\\left(f\\right)=\\frac{{SBSR}_{m,i}\\left(f\\right)/{SBSR}_{n,i}\\left(f\\right)}{{SBSR}_{m,j}\\left(f\\right)/{SBSR}_{n,j}\\left(f\\right)}$$\u003c/div\u003e\u003cdiv class=\"EquationNumber\"\u003e7\u003c/div\u003e\u003c/div\u003e\u003c/p\u003e\u003cp\u003eThis final expression reveals the physical meaning of the DTR: it compares the interevent SBSR fluctuation at station \u003cem\u003ei\u003c/em\u003e to that at station \u003cem\u003ej\u003c/em\u003e. If these fluctuations are driven by common path effects, the fluctuations at both stations should be highly correlated, and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:DTR\\left(f\\right)\\)\u003c/span\u003e\u003c/span\u003e should approach unity. Conversely, if the fluctuations are driven by site-specific factors (e.g., 2D/3D scattering effects dependent on the arrival direction), they will be largely independent between the two stations, and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:DTR\\left(f\\right)\\)\u003c/span\u003e\u003c/span\u003e will deviate significantly from unity.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec9\" class=\"Section2\"\u003e\u003ch2\u003e3.2. Analysis Procedure\u003c/h2\u003e\u003cp\u003eThe selection of station pairs for the final analysis was governed by two criteria designed to balance statistical robustness with physical constraints. First, the horizontal interstation distance had to be less than 20 km. This threshold represents a practical balance between securing a robust number of pairs for statistical analysis and limiting significant decorrelation of the incoming wavefield. Second, to ensure statistical stability, we required that each station pair have a minimum of 20 common earthquake pairs for each of the two analysis datasets described below. After applying these filtering steps, a final analysis set of 293 station pairs was selected (Fig.\u0026nbsp;4(a)).\u003c/p\u003e\u003cp\u003eTo systematically investigate the cause of SBSR variability, two distinct datasets of earthquake pairs were established for TSR analysis. The criteria for these datasets were designed to test specific hypotheses regarding the origin of variability. An example of the geographical distribution of earthquakes used to obtain these pairs relative to the MYGH04\u0026ndash;MYGH13 station pair is shown in Fig.\u0026nbsp;4(b) and 4(c).\u003c/p\u003e\u003cp\u003e\u003cb\u003eCollocated Earthquake Pairs (Dataset A)\u003c/b\u003e: This dataset comprises pairs of earthquakes that satisfy two conditions: a hypocenter separation smaller than 5 km and a M\u003csub\u003eJ\u003c/sub\u003e difference smaller than 0.2. This strategy was designed to minimize variations in both the propagation path and source size, thereby isolating SBSR variability caused by differences in the source radiation pattern or the complex site response to slight variations in the incoming wavefield. An example of this earthquake selection is shown for the MYGH04\u0026ndash;MYGH13 station pair in Fig.\u0026nbsp;4(b).\u003c/p\u003e\u003cp\u003e\u003cb\u003eCommon-Distance, Different-Path Pairs (Dataset B)\u003c/b\u003e: This dataset includes pairs of earthquakes that satisfy three conditions: a difference in hypocentral distance to the station pair smaller than 5 km, a M\u003csub\u003eJ\u003c/sub\u003e difference smaller than 0.2, and a back-azimuth separation of at least 45 degrees. This strategy was designed to compare events that travel over different paths but exhibit similar overall travel distances. These criteria for M\u003csub\u003eJ\u003c/sub\u003e and hypocentral distance ensure that the source term and the geometric attenuation are comparable for both events in a pair. However, because the paths encompass different crustal volumes, the anelastic attenuation (Q value) is expected to differ. This dataset facilitates a direct assessment of whether the TSR method consistently captures these differing path effects at both the surface and the borehole, serving as validation of the overall approach. An example of this earthquake selection is shown for the MYGH04\u0026ndash;MYGH13 station pair in Fig.\u0026nbsp;4(c).\u003c/p\u003e\u003cp\u003eFor each of the 293 station pairs and for each of the two earthquake pairing strategies, two correlation analyses were conducted across the 1\u0026ndash;10 Hz frequency range. These analyses provided the quantitative metrics for this study.\u003c/p\u003e\u003cp\u003e\u003cstrong\u003ePath-Effect Correlation\u003c/strong\u003e\u003cp\u003ePearson\u0026rsquo;s correlation coefficient between the natural logarithm of the surface-TSR (ln(\u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{TSR}^{S}\\)\u003c/span\u003e\u003c/span\u003e)) and the borehole-TSR (ln(\u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{TSR}^{B}\\)\u003c/span\u003e\u003c/span\u003e​)) was calculated. This calculation serves as a validation or quality control step for the methodology. A high correlation coefficient indicates that the TSR method can effectively isolate the same path-related term at both the surface and the borehole.\u003c/p\u003e\u003c/p\u003e\u003cp\u003e\u003cstrong\u003eSite-Specific Variability Correlations\u003c/strong\u003e\u003cp\u003ePearson\u0026rsquo;s correlation coefficient between the logarithmic interevent SBSR fluctuations at the two sites (between ln[\u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{SBSR}_{m,i}\\left(f\\right)/{SBSR}_{n,i}\\left(f\\right)\\)\u003c/span\u003e\u003c/span\u003e] and ln[\u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{SBSR}_{m,j}\\left(f\\right)/{SBSR}_{n,j}\\left(f\\right)\\)\u003c/span\u003e\u003c/span\u003e​]) was calculated. This represents the primary metric of the study. It directly reflects the degree of spatial coherence in the observed site response variability. A low correlation coefficient suggests that the event-to-event fluctuations in the SBSR at one station are independent of those at a nearby station, strongly suggesting that variability is controlled by local, site-specific conditions.\u003c/p\u003e\u003c/p\u003e\u003c/div\u003e"},{"header":"4. Results and Discussion","content":"\u003cp\u003eThis section presents the results of the modified TSR approach, quantitatively assessing the spatial coherence of event-to-event SBSR variability. The analysis is presented first for the overall statistical distributions of all 293 station pairs, followed by a detailed examination of a representative station pair (MYGH04\u0026ndash;MYGH13) to illustrate the underlying physical mechanisms.\u003c/p\u003e\u003cdiv id=\"Sec11\" class=\"Section2\"\u003e\u003ch2\u003e4.1. Path-Effect Correlation: Surface-TSR versus Borehole-TSR\u003c/h2\u003e\u003cp\u003eThe correlation between the surface-TSR (ln(\u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{TSR}^{S}\\)\u003c/span\u003e\u003c/span\u003e)) and the borehole-TSR (ln(\u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{TSR}^{B}\\)\u003c/span\u003e\u003c/span\u003e)) serves as a critical validation metric for our approach. A high correlation indicates that the TSR successfully isolates the same differential path-effect term at both the surface and borehole sensors, regardless of the complex wavefield in the borehole. The frequency-dependent Pearson correlation coefficients for this path-effect term are shown in Fig.\u0026nbsp;5. The results for the collocated earthquake pairs (Dataset A) are provided in the lower panel of Fig.\u0026nbsp;5(a), while the results for the common-distance, different-path pairs (Dataset B) are shown in the lower panel of Fig.\u0026nbsp;5(b).\u003c/p\u003e\u003cp\u003eFor both datasets, the overall statistical distributions demonstrate strong positive correlations between the surface-TSR and the borehole-TSR. In the case of the collocated pairs in Fig.\u0026nbsp;5(a) (lower panel), the results of the analysis show that the median correlation coefficient is approximately 0.9 at 1 Hz, and while it gradually decreases with increasing frequency, it remains above 0.5 up to 10 Hz. This finding indicates that for more than half of the station pairs, the path-related spectral modulations are captured consistently at both sensor depths. The results for the different-path dataset in Fig.\u0026nbsp;5(b) (lower panel) reveal a very similar trend, with a robust median correlation across the entire frequency band.\u003c/p\u003e\u003cp\u003eThis high degree of correlation across a large and diverse set of station pairs confirms the fundamental premise of the analysis: the path-effect term, \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{P}_{m,i}\\left(f\\right)/{P}_{n,i}\\left(f\\right)\\)\u003c/span\u003e\u003c/span\u003e, is largely independent of the local site effect and is consistently captured at both the surface and the borehole. This consistency holds even when the paths substantially differ, as in Dataset B, which suggests that the anelastic attenuation characteristics (related to the Q value) along different paths are robustly captured by the TSR. The successful isolation of this common path term is essential, as it allows us to attribute the remaining variability to other factors.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec12\" class=\"Section2\"\u003e\u003ch2\u003e4.2. Site-Specific Variability Correlation: Interstation SBSR Fluctuations\u003c/h2\u003e\u003cp\u003eHaving validated the method's ability to isolate a common path term in Section \u003cspan refid=\"Sec11\" class=\"InternalRef\"\u003e4.1\u003c/span\u003e, the core of our subsequent analysis lies in determining whether the event-to-event SBSR variability is a spatially coherent effect driven by regional path phenomena, or a spatially incoherent effect driven by local site conditions. We therefore calculated the correlation of the interevent SBSR fluctuations (\u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{SBSR}_{m,i}\\left(f\\right)/{SBSR}_{n,i}\\left(f\\right)\\)\u003c/span\u003e\u003c/span\u003e) between each of the 293 station pairs. The results are shown in the upper panels of Fig.\u0026nbsp;5(a) and 5(b) and are consistent across both datasets. The median interstation correlation of the SBSR fluctuations remains close to zero across the entire 1\u0026ndash;10 Hz frequency range. The interquartile range (IQR), indicated by the dark shaded area, is also narrowly centered on zero, typically ranging from approximately \u0026minus;\u0026thinsp;0.3 to +\u0026thinsp;0.3. This near-zero correlation indicates a lack of spatial coherence in SBSR variability. If the azimuthal dependence of the SBSR was caused primarily by large-scale, path-dependent effects such as directional variations in regional attenuation or scattering from deep crustal heterogeneities, one would expect to observe a positive correlation between nearby stations, as they would share a significant portion of their propagation paths.\u003c/p\u003e\u003cp\u003eThe absence of such a correlation strongly suggests that the observed event-to-event fluctuations are controlled primarily by factors unique to each individual site. These site-specific complexities, such as small-scale geological irregularities, surface topography, or basin-edge effects, are known to interact with the incoming wavefield in a manner that is highly sensitive to the arrival direction (Bard and Bouchon, \u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e1985\u003c/span\u003e). The resulting scattered wavefield and directional resonance patterns are therefore unique to the immediate vicinity of each station, leading to the observed spatial decorrelation of SBSR variability even between stations separated by less than 20 km.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec13\" class=\"Section2\"\u003e\u003ch2\u003e4.3. Case Study: Influence of Azimuthal Dependence at MYGH04\u0026ndash;MYGH13\u003c/h2\u003e\u003cp\u003eThe overall statistical results in Section \u003cspan refid=\"Sec12\" class=\"InternalRef\"\u003e4.2\u003c/span\u003e indicate that SBSR variability is a spatially incoherent phenomenon controlled by local site conditions. To explore the physical mechanisms underlying this general conclusion in greater detail, we examined the results for the MYGH04\u0026ndash;MYGH13 station pair, highlighted by magenta lines in Fig.\u0026nbsp;5. As established in the preliminary analysis (Fig.\u0026nbsp;2(b)), the MYGH13 station exhibits exceptionally high azimuthal dependence of its SBSR.\u003c/p\u003e\u003cp\u003eLet us first consider the results for Dataset B (Fig.\u0026nbsp;5(b)), which represents the more complex scenario. The path-effect correlation for this pair (lower panel, magenta line) deviates significantly from the median trend. It indicates notable decreases in correlation, decreasing to nearly zero within the 3\u0026ndash;4.5 Hz and 6\u0026ndash;9.5 Hz frequency bands. These frequency bands coincide precisely with the frequencies at which the SBSR at MYGH13 exhibits its largest peak shifts and shape changes with back-azimuth. This result is instructive: it demonstrates that when site effects highly depend on the arrival direction, they are no longer completely canceled during TSR calculation. The site term becomes event-dependent (i.e., \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{G}_{m,i}^{S}\\left(f\\right)\\ne\\:{G}_{n,i}^{S}\\left(f\\right)\\)\u003c/span\u003e\u003c/span\u003e or \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{G}_{m,i}^{B}\\left(f\\right)\\ne\\:{G}_{n,i}^{B}\\left(f\\right)\\)\u003c/span\u003e\u003c/span\u003e), violating the assumption underlying the standard TSR method.\u003c/p\u003e\u003cp\u003eThe violation means the site-specific response is not fully canceled and instead distorts the calculated path-effect term, which reduces the correlation between the surface-TSR and the borehole-TSR, directly influencing the interstation correlation of SBSR fluctuations (upper panel, magenta line). Within the same frequency bands (3\u0026ndash;4 Hz and 6\u0026ndash;9.5 Hz) where the path-effect correlation is low, the SBSR fluctuation correlation for MYGH04\u0026ndash;MYGH13 becomes strongly negative, approaching \u0026minus;\u0026thinsp;0.5. This finding indicates an apparent inverse correlation between the SBSR fluctuations at the two sites. However, this inverse correlation should be interpreted with caution. It occurs precisely where the underlying assumption of the analysis, namely, that the common path effect is successfully isolated and removed, breaks down.\u003c/p\u003e\u003cp\u003eGiven that this negative correlation occurs precisely where the underlying assumption of the analysis\u0026mdash;that the common path effect is successfully isolated\u0026mdash;breaks down, we interpret this result as a methodological artifact rather than a true physical phenomenon. The inclusion of such strong negative correlations from a subset of station pairs explains the greater overall statistical dispersion observed in Dataset B, which is reflected by the wider whisker range in the upper panel of Fig.\u0026nbsp;5(b) compared with that of Fig.\u0026nbsp;5(a). Within the frequency bands where the path-effect correlation is high (i.e., below 3 Hz and between 4.5 and 6 Hz), the interstation SBSR fluctuation correlation is indeed close to zero.\u003c/p\u003e\u003cp\u003eThe results for Dataset A (Fig.\u0026nbsp;5(a)) provide a clearer picture than do those for Dataset B. For Dataset A, the path-effect correlation for MYGH04\u0026ndash;MYGH13 (lower panel, magenta line) is consistently high, following the median trend. More importantly, the interstation correlation of SBSR fluctuations (upper panel, magenta line) remains consistently close to zero across the entire frequency range. Via the use of earthquake pairs where path-related differences are nearly eliminated, we isolate the site response to slight variations in the incoming wavefield. The fact that the resulting SBSR fluctuations at station MYGH13 are not correlated with those at the nearby MYGH04 station, even in frequency bands susceptible to notable azimuthal effects, demonstrates that this variability is a local site-specific effect.\u003c/p\u003e\u003cp\u003eThe analysis of the MYGH04\u0026ndash;MYGH13 station pair provides two important observations. First, the negative correlation of SBSR fluctuations for Dataset B coincides with the frequency bands where the path-effect correlation is low. This suggests that the negative correlation is an artifact of applying the TSR method, in which the underlying assumptions are violated by notable, azimuth-dependent site effects. Second, the analysis of Dataset A, where path effects are minimal, reveals a near-zero correlation of SBSR fluctuations across all frequencies. These results indicate that the considerable azimuthal variability at MYGH13 is a localized characteristic attributable to the site conditions and is not spatially correlated with the nearby station.\u003c/p\u003e\u003c/div\u003e"},{"header":"5. Conclusion","content":"\u003cp\u003eThis study investigated the significant, azimuth-dependent variability of SBSRs using a novel TSR-based approach to distinguish between path- specific and site-specific effects. Our analysis of 293 KiK-net station pairs revealed a high correlation between path-effect terms derived from surface and borehole records but a near-zero interstation correlation of event-to-event SBSR fluctuations.\u003c/p\u003e\u003cp\u003eThe observed spatial decorrelation of SBSR variability provides strong evidence that its origin is predominantly local and site-specific. If this variability were driven by regional-scale path effects, a significant positive correlation would be expected between nearby stations. The absence of such a correlation indicates that the azimuthal dependence of site response is governed by site-specific complexities, such as 2D/3D geological structures.\u003c/p\u003e\u003cp\u003eFurthermore, our experimental design, utilizing two distinct datasets, offers a potential framework for decomposing site-specific uncertainty. The variability observed using collocated earthquake pairs (Dataset A) can be interpreted as representing the aleatory uncertainty, the inherent randomness in the site response. In contrast, the variability from the different-azimuth pairs (Dataset B) encompasses this aleatory uncertainty plus the epistemic uncertainty associated with the systematic, unmodeled effects of azimuthal dependence. This suggests a pathway to validate the appropriateness and potential impact of applying non-ergodic site terms, by demonstrating how much of the total site-specific variability they can successfully explain.\u003c/p\u003e\u003cp\u003eThese results have direct implications for the next generation of seismic hazard assessment. By reclassifying a significant component of ground motion variability from random uncertainty to a deterministic, local site attribute, this work provides the physical basis and empirical data needed to develop more robust non-ergodic ground motion models (GMMs). Modern GMMs increasingly seek to break the ergodic assumption by separating the total variability into regional and site-specific components (Shahi and Baker, \u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e2014\u003c/span\u003e; Baker and Bradley, 2014). This finding directly supports the ongoing development of non-ergodic GMMs. While these models already aim to isolate site-specific variability (Shahi and Baker, \u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e2014\u003c/span\u003e), the indeterminate origin of its azimuthal component has been a key source of epistemic uncertainty. By demonstrating that this azimuthal dependence is a stable, site-specific phenomenon, our work provides a physical validation for treating this variability as a deterministic local attribute, paving the way for its explicit incorporation into future models.\u003c/p\u003e\u003cp\u003eWhile this study provides a new framework for understanding site-specific variability, it is important to acknowledge its limitations. First, our methodology, while effective at identifying the origin of azimuthal dependence (site versus path), does not fully characterize the specific physical mechanisms (e.g., scattering from topography versus basin resonance) at each individual site. Second, our findings are based on the dense KiK-net in Japan, and their generalizability to other tectonic regions requires further investigation. Furthermore, our findings are based on the dense KiK-net in Japan; their generalizability to other tectonic regions requires further investigation. Future work should therefore focus on applying this methodology to other regions and integrating the results with detailed geological and topographical data to build more quantitative models of these 2D/3D effects. Such efforts are essential for developing the practical parameterizations needed to refine non-ergodic models (Lavrentiadis et al., \u003cspan citationid=\"CR10\" class=\"CitationRef\"\u003e2023\u003c/span\u003e) and will ultimately lead to more reliable PSHA.\u003c/p\u003e"},{"header":"Abbreviations","content":"\u003cdiv class=\"DefinitionList\"\u003e\u003cdiv class=\"DefinitionListEntry\"\u003e\u003cdiv class=\"Term\"\u003e1D/2D/3D\u003c/div\u003e\u003cdiv class=\"Description\"\u003e\u003cp\u003eOne-/two-/three-dimensional\u003c/p\u003e\u003c/div\u003e\u003c/div\u003e\u003cdiv class=\"DefinitionListEntry\"\u003e\u003cdiv class=\"Term\"\u003eCSMIP\u003c/div\u003e\u003cdiv class=\"Description\"\u003e\u003cp\u003eCalifornia Strong-Motion Instrumentation Program\u003c/p\u003e\u003c/div\u003e\u003c/div\u003e\u003cdiv class=\"DefinitionListEntry\"\u003e\u003cdiv class=\"Term\"\u003eDTR\u003c/div\u003e\u003cdiv class=\"Description\"\u003e\u003cp\u003eDual twofold spectral ratio\u003c/p\u003e\u003c/div\u003e\u003c/div\u003e\u003cdiv class=\"DefinitionListEntry\"\u003e\u003cdiv class=\"Term\"\u003eETF\u003c/div\u003e\u003cdiv class=\"Description\"\u003e\u003cp\u003eEmpirical transfer function\u003c/p\u003e\u003c/div\u003e\u003c/div\u003e\u003cdiv class=\"DefinitionListEntry\"\u003e\u003cdiv class=\"Term\"\u003e\u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:F\\)\u003c/span\u003e\u003c/span\u003e\u003c/div\u003e\u003cdiv class=\"Description\"\u003e\u003cp\u003eFourier spectra of the observation records\u003c/p\u003e\u003c/div\u003e\u003c/div\u003e\u003cdiv class=\"DefinitionListEntry\"\u003e\u003cdiv class=\"Term\"\u003e\u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:G\\)\u003c/span\u003e\u003c/span\u003e\u003c/div\u003e\u003cdiv class=\"Description\"\u003e\u003cp\u003eGround amplification function\u003c/p\u003e\u003c/div\u003e\u003c/div\u003e\u003cdiv class=\"DefinitionListEntry\"\u003e\u003cdiv class=\"Term\"\u003eGMM\u003c/div\u003e\u003cdiv class=\"Description\"\u003e\u003cp\u003eGround motion model\u003c/p\u003e\u003c/div\u003e\u003c/div\u003e\u003cdiv class=\"DefinitionListEntry\"\u003e\u003cdiv class=\"Term\"\u003eIQR\u003c/div\u003e\u003cdiv class=\"Description\"\u003e\u003cp\u003eInterquartile range\u003c/p\u003e\u003c/div\u003e\u003c/div\u003e\u003cdiv class=\"DefinitionListEntry\"\u003e\u003cdiv class=\"Term\"\u003eJMA\u003c/div\u003e\u003cdiv class=\"Description\"\u003e\u003cp\u003eJapan Meteorological Agency\u003c/p\u003e\u003c/div\u003e\u003c/div\u003e\u003cdiv class=\"DefinitionListEntry\"\u003e\u003cdiv class=\"Term\"\u003eKiK-net\u003c/div\u003e\u003cdiv class=\"Description\"\u003e\u003cp\u003eKiban\u0026ndash;Kyoshin network\u003c/p\u003e\u003c/div\u003e\u003c/div\u003e\u003cdiv class=\"DefinitionListEntry\"\u003e\u003cdiv class=\"Term\"\u003eln\u003c/div\u003e\u003cdiv class=\"Description\"\u003e\u003cp\u003eNatural logarithm\u003c/p\u003e\u003c/div\u003e\u003c/div\u003e\u003cdiv class=\"DefinitionListEntry\"\u003e\u003cdiv class=\"Term\"\u003eMANOVA\u003c/div\u003e\u003cdiv class=\"Description\"\u003e\u003cp\u003eMultivariate analysis of variance\u003c/p\u003e\u003c/div\u003e\u003c/div\u003e\u003cdiv class=\"DefinitionListEntry\"\u003e\u003cdiv class=\"Term\"\u003eM\u003csub\u003eJ\u003c/sub\u003e\u003c/div\u003e\u003cdiv class=\"Description\"\u003e\u003cp\u003eMagnitude defined by the JMA\u003c/p\u003e\u003c/div\u003e\u003c/div\u003e\u003cdiv class=\"DefinitionListEntry\"\u003e\u003cdiv class=\"Term\"\u003eNIED\u003c/div\u003e\u003cdiv class=\"Description\"\u003e\u003cp\u003eNational Research Institute for Earth Science and Disaster Resilience\u003c/p\u003e\u003c/div\u003e\u003c/div\u003e\u003cdiv class=\"DefinitionListEntry\"\u003e\u003cdiv class=\"Term\"\u003e\u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:P\\)\u003c/span\u003e\u003c/span\u003e\u003c/div\u003e\u003cdiv class=\"Description\"\u003e\u003cp\u003ePath transfer function\u003c/p\u003e\u003c/div\u003e\u003c/div\u003e\u003cdiv class=\"DefinitionListEntry\"\u003e\u003cdiv class=\"Term\"\u003ePGA\u003c/div\u003e\u003cdiv class=\"Description\"\u003e\u003cp\u003ePeak ground acceleration\u003c/p\u003e\u003c/div\u003e\u003c/div\u003e\u003cdiv class=\"DefinitionListEntry\"\u003e\u003cdiv class=\"Term\"\u003ePSHA\u003c/div\u003e\u003cdiv class=\"Description\"\u003e\u003cp\u003eProbabilistic seismic hazard analysis\u003c/p\u003e\u003c/div\u003e\u003c/div\u003e\u003cdiv class=\"DefinitionListEntry\"\u003e\u003cdiv class=\"Term\"\u003eRSF\u003c/div\u003e\u003cdiv class=\"Description\"\u003e\u003cp\u003eRelative site factor\u003c/p\u003e\u003c/div\u003e\u003c/div\u003e\u003cdiv class=\"DefinitionListEntry\"\u003e\u003cdiv class=\"Term\"\u003e\u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:S\\)\u003c/span\u003e\u003c/span\u003e\u003c/div\u003e\u003cdiv class=\"Description\"\u003e\u003cp\u003eSource spectra\u003c/p\u003e\u003c/div\u003e\u003c/div\u003e\u003cdiv class=\"DefinitionListEntry\"\u003e\u003cdiv class=\"Term\"\u003eSBSR\u003c/div\u003e\u003cdiv class=\"Description\"\u003e\u003cp\u003eSurface-to-borehole spectral ratio\u003c/p\u003e\u003c/div\u003e\u003c/div\u003e\u003cdiv class=\"DefinitionListEntry\"\u003e\u003cdiv class=\"Term\"\u003eSD\u003c/div\u003e\u003cdiv class=\"Description\"\u003e\u003cp\u003eStandard deviation\u003c/p\u003e\u003c/div\u003e\u003c/div\u003e\u003cdiv class=\"DefinitionListEntry\"\u003e\u003cdiv class=\"Term\"\u003eTSR\u003c/div\u003e\u003cdiv class=\"Description\"\u003e\u003cp\u003eTwofold spectral ratio\u003c/p\u003e\u003c/div\u003e\u003c/div\u003e\u003cdiv class=\"DefinitionListEntry\"\u003e\u003cdiv class=\"Term\"\u003eTTF\u003c/div\u003e\u003cdiv class=\"Description\"\u003e\u003cp\u003eTheoretical transfer function\u003c/p\u003e\u003c/div\u003e\u003c/div\u003e\u003c/div\u003e"},{"header":"Declarations","content":"\u003cp\u003e\u003cstrong\u003eEthics approval and consent to participate\u003c/strong\u003e\u003cp\u003eNot applicable\u003c/p\u003e\u003cp\u003e\u003cstrong\u003eConsent for publication\u003c/strong\u003e\u003cp\u003eNot applicable\u003c/p\u003e\u003c/p\u003e\u003ch2\u003eAuthors' contributions\u003c/h2\u003e\u003cp\u003eThe author conceived the study, analyzed the data, and wrote the manuscript. Dr. Yamamoto, Dr. Tokumitsu, and Dr. Uchiyama contributed to the development of the research direction through regular discussions and provided critical feedback on the manuscript. Prof. Takai, as the research supervisor, provided overall guidance for this research and reviewed the manuscript. All the authors have read and approved the final manuscript.\u003c/p\u003e\u003ch2\u003eAcknowledgements\u003c/h2\u003e\u003cp\u003eWe thank the NIED for providing the KiK-net observation records.\u003c/p\u003e\u003ch2\u003eAvailability of data and materials\u003c/h2\u003e\u003cp\u003eThe strong-motion data used in this study were obtained from the Kiban\u0026ndash;Kyoshin network (KiK-net), operated by the National Research Institute for Earth Science and Disaster Resilience (NIED). The data are publicly available from the NIED website at \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://www.kyoshin.bosai.go.jp/\u003c/span\u003e\u003cspan address=\"https://www.kyoshin.bosai.go.jp/\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e. The analysis code used in this study is available from the corresponding author upon reasonable request.\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\u003cli\u003e\u003cspan\u003eAki K, Richards PG (2002) Quantitative Seismology, 2nd edn. University Science Books\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eArchuleta RJ, Seale SH, Sangas PV, Baker LM, Swain ST (1992) Garner Valley downhole array of accelerometers: Instrumentation and preliminary data analysis. Bull Seismol Soc Am 82:1592\u0026ndash;1621. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1785/BSSA0820041592\u003c/span\u003e\u003cspan address=\"10.1785/BSSA0820041592\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eBaker JW, Bradley BA (2017) Intensity Measure Correlations Observed in the NGA-West2 Database, and Dependence of Correlations on Rupture and Site Parameters. Earthq Spectra 33(1):145\u0026ndash;156. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1193/060716eqs095m\u003c/span\u003e\u003cspan address=\"10.1193/060716eqs095m\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eBard PY, Bouchon M (1985) The two-dimensional resonance of sediment-filled valleys. Bull Seismol Soc Am 75(2):519\u0026ndash;541. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1785/BSSA0750020519\u003c/span\u003e\u003cspan address=\"10.1785/BSSA0750020519\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eField EH, Jacob KH (1995) A comparison and test of various site-response estimation techniques, including three that are not reference-site dependent. Bull Seismol Soc Am 85:1127\u0026ndash;1143. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1785/BSSA0850041127\u003c/span\u003e\u003cspan address=\"10.1785/BSSA0850041127\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eHair JF Jr, Black WC, Babin B, Anderson RE (2018) Multivariate Data Analysis, Eighth Edition. Cengage Learning EMEA\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eIkeura T (2020) A method for investigating the attenuation characteristics of strong motions without attenuation functions: Application to strong motions in the Tohoku and Kanto regions caused by M7 class plate boundary earthquakes. J Japan Association Earthq Eng. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.5610/jaee.20.7_133\u003c/span\u003e\u003cspan address=\"10.5610/jaee.20.7_133\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e(in Japanese) 20:7_133-7_157\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eIzumitani Y (2000) Qs-value in southern Kyushu evaluated from double spectral ratio of strong motion records. Journal of Japan Society of Civil Engineers, 2000:225\u0026ndash;230. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.2208/jscej.2000.640_225\u003c/span\u003e\u003cspan address=\"10.2208/jscej.2000.640_225\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (in Japanese)\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eJoyner WB, Warrick RE, Oliver AA III (1976) Analysis of seismograms from a downhole array in sediments near San Francisco Bay. Bull Seismol Soc Am 66(3):937\u0026ndash;958. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1785/BSSA0660030937\u003c/span\u003e\u003cspan address=\"10.1785/BSSA0660030937\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eLavrentiadis G, Abrahamson NA, Kuehn NM (2023) A non-ergodic effective amplitude ground-motion model for California. Bull Earthq Eng 21:5233\u0026ndash;5264. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1007/s10518-021-01206-w\u003c/span\u003e\u003cspan address=\"10.1007/s10518-021-01206-w\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eMatsuzawa T, Takeo M, Ide S, Iio Y, Ito H, Imanishi K, Horiuchi S (2004) Estimation of the S-wave attenuation in the western Nagano region from twofold spectral ratio. J Seismological Soc Japan Ser 2 56:75\u0026ndash;88. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.4294/zisin1948.56.1_75\u003c/span\u003e\u003cspan address=\"10.4294/zisin1948.56.1_75\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e(in Japanese)\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eNational Research Institute for Earth Science and Disaster Resilience (2019) NIED K-NET, KiK-net. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.17598/NIED.0004\u003c/span\u003e\u003cspan address=\"10.17598/NIED.0004\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eNoguchi S, Sasatani T (2011) Nonlinear soil response and its effects on strong ground motions during the 2003 Miyagi-Oki intraslab earthquake. J Seismological Soc Japan Ser 2 63:165\u0026ndash;187. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.4294/zisin.63.165\u003c/span\u003e\u003cspan address=\"10.4294/zisin.63.165\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e(in Japanese)\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eParker GA, Stewart JP, Boore DM, Atkinson GM, Hassani B (2021) NGA-subduction global ground motion models with regional adjustment factors. Earthq Spectra 38:456\u0026ndash;493. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1177/87552930211034889\u003c/span\u003e\u003cspan address=\"10.1177/87552930211034889\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003ePedregosa F, Varoquaux G, Gramfort A, Michel V, Thirion B, Grisel O, Blondel M, Prettenhofer P, Weiss R, Dubourg V, Vanderplas J, Passos A, Cournapeau D, Brucher M, Perrot M, Duchesnay \u0026Eacute; (2011) Scikit-learn: Machine Learning in Python. J Mach Learn Res 12:2825\u0026ndash;2830\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003ePilz M, Cotton F (2019) Does the one-dimensional assumption hold for site response analysis? A study of seismic site responses and implication for ground motion assessment using KiK-Net strong-motion data. Earthq Spectra 35:883\u0026ndash;905. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1193/050718EQS113M\u003c/span\u003e\u003cspan address=\"10.1193/050718EQS113M\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eSeabold S, Perktold J (2010) statsmodels: Econometric and statistical modeling with python. In: Proceedings of the 9th Python in Science Conference\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eSgobba S, Felicetta C, Bortolotti T, Menafoglio A, Lanzano G, Pacor F (2024) A geostatistical modelling of empirical amplification functions and related site proxies for shaking scenarios in central Italy. Soil Dyn Earthq Eng 179:108496. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1016/j.soildyn.2024.108496\u003c/span\u003e\u003cspan address=\"10.1016/j.soildyn.2024.108496\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eShahi SK, Baker JW (2014) NGA-West2 Models for Ground Motion Directionality. Earthq Spectra 30(3):1285\u0026ndash;1300. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1193/040913EQS097M\u003c/span\u003e\u003cspan address=\"10.1193/040913EQS097M\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eShiomi K (2012) An improved method for estimating the installation azimuth information of NIED Hi-net borehole seismometers. Rep Natl Res Inst Earth Sci Disaster Resil 80:1\u0026ndash;20 (in Japanese)\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eShima E (1962) Modification of Seismic Waves in Superficial Layers as Verified by Comparative Observations on and Beneath the Surface. Bull Earthq Res Inst Tokyo 40:187\u0026ndash;259\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eTao Y, Rathje E (2020) The importance of distinguishing pseudoresonances and outcrop resonances in downhole array data. Bull Seismol Soc Am 110:288\u0026ndash;294. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1785/0120190097\u003c/span\u003e\u003cspan address=\"10.1785/0120190097\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eThompson EM, Baise LG, Tanaka Y, Kayen RE (2012) A taxonomy of site response complexity. Soil Dyn Earthq Eng 41:32\u0026ndash;43. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1016/j.soildyn.2012.04.005\u003c/span\u003e\u003cspan address=\"10.1016/j.soildyn.2012.04.005\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eUeno H, Hatakeyama S, Aketagawa T, Funazaki A, Hamada N (2002) Improvement of hypocenter determination procedures in the Japan Meteorological Agency. Seismological Bull Japan Meteorological Agency 65:123\u0026ndash;134 (in Japanese)\u003c/span\u003e\u003c/li\u003e\u003c/ol\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":false,"hideJournal":true,"highlight":"","institution":"","isAcceptedByJournal":false,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true},"keywords":"Surface-to-Borehole Spectral Ratio, Azimuthal Variation, Twofold Spectral Ratio, Site Amplification, Site-Specific Variability, KiK-net","lastPublishedDoi":"10.21203/rs.3.rs-7806012/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-7806012/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003cp\u003eIn this study, we show that the widely observed azimuthal dependence of the surface-to-borehole spectral ratios (SBSRs) is governed primarily by local, site-specific complexities rather than regional path-specific effects, resolving a fundamental ambiguity in site response analysis. This ambiguity arises from the significant event-to-event variability seen in empirical data, which contradicts the stable, time-invariant theoretical transfer function predicted by one-dimensional (1D) site response analysis and introduces considerable epistemic uncertainty into ground motion models (GMMs). A key challenge has been the indeterminate origin of this variance whether attributable to the propagation path or the local site. To systematically decompose these contributing effects, we developed an analytical framework based on the twofold spectral ratio (TSR) method. Applying this framework to a large dataset from 293 nearby vertical array station pairs from Japan's Kiban\u0026ndash;Kyoshin network (KiK-net), our analysis was designed to isolate and nullify the common path-effect term, enabling a direct, quantitative evaluation of site-specific variability. The analysis yielded two empirical results. First, the path-effect term isolated by the TSR showed high correlation between surface and borehole records, with the median correlation coefficient remaining above 0.5 across the entire 1\u0026ndash;10 Hz frequency band, validating the successful separation of this component. Second, the residual event-to-event SBSR fluctuations between adjacent stations exhibited a near-zero spatial correlation across the entire frequency range, which is inconsistent with the hypothesis of a shared, regional cause. This reclassification of variability from a random uncertainty to a deterministic, site-specific attribute has direct implications for advancing seismic hazard assessment. Our dual-dataset approach provides a framework to validate the decomposition of site-specific variability into its aleatory and epistemic components. This provides the physical basis for developing non-ergodic GMMs, which can reduce overall uncertainty by explicitly modeling the predictable epistemic component, facilitating more sophisticated, physics-based site characterization.\u003c/p\u003e","manuscriptTitle":"Investigating the Origin of Azimuthal Site Response Variability using KiK-net Data: A Twofold Spectral Ratio Approach","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2025-11-07 06:23:45","doi":"10.21203/rs.3.rs-7806012/v1","editorialEvents":[{"type":"communityComments","content":0}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true}}],"origin":"","ownerIdentity":"ea255ea0-4850-4988-a82d-a766e6977590","owner":[],"postedDate":"November 7th, 2025","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"posted","subjectAreas":[],"tags":[],"updatedAt":"2026-05-18T08:21:50+00:00","versionOfRecord":[],"versionCreatedAt":"2025-11-07 06:23:45","video":"","vorDoi":"","vorDoiUrl":"","workflowStages":[]},"version":"v1","identity":"rs-7806012","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-7806012","identity":"rs-7806012","version":["v1"]},"buildId":"8U1c8b4HqxoKbykW_rLl7","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}