Horizontal seismic wave at ground surface from transfer function based on ambient noise

Earthquake detection can be improved by ensuring that seismometer sites experience little artificial noise in the surrounding environment. To minimize noise, seismological stations should be positioned in rocky mountainous areas without nearby valleys, away from significant human activity. However, such surface sites may be scarce when constructing dense monitoring networks, necessitating the use of underground sites to ensure low noise levels. The Korean Meteorological Administration is currently installing new underground seismometers to increase seismic monitoring capacity. However, seismic data on the ground surface are also required for engineering technological developments (to reduce damage to structural components). Therefore, borehole seismic stations without surface seismometers need to estimate ground surface motion from borehole record data. We propose a transfer function that converts motion within boreholes to surface seismic waves using ambient noise, thereby facilitating estimation of ground surface motions using borehole seismometers. As a result, predicting ground surface motion from borehole record data becomes possible.


Introduction
The Korea Meteorological Administration (KMA) detects earthquakes and analyzes the source information using its seismic monitoring network. The Korean seismic monitoring network began with two stations in 1978, which underwent rapid modernization in the 1990s (KMA, 2001;Park et al., 2011;Cho et al., 2022). The number of stations further increased following the 2016 Gyeongju earthquake (Cho et al., 2022). Currently, 282 seismometers are installed, with borehole-type seismometers being applied to new or renovated stations since 2016. Recently, seismometers, have been installed in subsurface bedrock layers, to reduce the influence of day-to-day noise. Due to their low ambient noise, these borehole stations can facilitate outstanding detection of micro earthquakes and the first signs of P-waves. As of May 2022, borehole and surface seismometers account for 85% and 15%, respectively, of seismic observation stations.
With regard to earthquake early warning (EEW) systems, installation of borehole stations is highly efficient (Cho et al., 2022;Jang et al., 2023). However, as we aimed to assess earthquake-induced damage from ground motion, we cannot neglect the management of recordings on the ground surface. Unfortunately, high installation costs and budgetary limitations have restricted the establishment of both surface and borehole observatories. As a result, there are only limited studies of how vibrations are transferred from boreholes to the ground surface in the context of the topography of the Korean Peninsula.
To solve this problem, a few studies have applied a site amplitude coefficient (SAC) (Borcherdt, 1994;BSSC, 1997;Dobry et al., 2000) that has been used in the western regions of the United States. However, this method may not be appropriate, because the geological and topographical characteristics of Korea and the United States differ, and therefore a detailed review is needed to apply this method to the Korean Peninsula (Manandhar et al., 2018). In addition, estimation techniques related to earthquake damage are based on an empirical formula that utilizes the amplitude of ground surface vibrations; thus, it is crucial to determine surface vibrations (McCann et al., 1980;Shinozuka et al., 2000;Padgett and DesRoches, 2008). To this end, it would be highly effective to install surface sensors in the borehole stations; however, this was not possible because of the available budget in Korea. Thus, it is necessary to develop a technique to estimate ground surface vibrations using subsurface sensors in order to determine ground surface behavior.
Two approaches can be used to estimate ground surface vibrations. First, in the SAC method, which was first introduced in seismic design, the underlying hard bedrock conditions are assumed to match those of rock outcrops. Based on this, an increment-decrement coefficient is applied according to the average shear-wave velocity up to 30 m from the target ground. The modified Mercalli intensity (MMI) system of the KMA uses the SAC of the Borcherdt model (Borcherdt, 1994). However, the site coefficient is not suitable as a ground surface prediction model for borehole stations, as the within-motion in the rock around the borehole differs from the rock outcrop motion in the free field (Kwak et al., 2022;Lai et al., 2022). Second, the ground response can be interpreted via numerical analysis. This approach requires characterization of materials comprising the thin soil layer. This can be achieved through ground surveys such as suspension PS (S-PS) logging, downhole testing, multichannel analysis of surface waves (MASW), and standard penetration tests. This ground surface analysis method has been shown to be highly effective in predicting ground surface vibrations (Park and Hashash, 2004). Kwak et al. (2022) proposed estimation of the ground surface response by combining the ground motion model (GMM) and 1D site response. However, this required a prior site characterize. It can be prohibitively expensive to perform such surveys at a large number of stations; moreover, at some sites, it is not feasible to conduct surveys at all. Therefore, a new estimation method without site characterization is required to predict ground surface vibrations at borehole stations. We tried to use ambient noise for the non-experimental approaches. At this time, an ambient noise signal is utilized because it includes a resonance effect on the media. Many studies have already analyzed the correlation with the predominant and/or fundamental frequency of the site characterization from ambient noise (Albarello and Lunedei, 2010;García-Fernández and Jiménez, 2012;Cultrera et al., 2014;Moisidi et al., 2015;Asten and Hayashi, 2018). The aim of this study was to determine a transfer function to convert motion recorded at a borehole to that occurring at the ground surface, without the need for site-specific response analysis. To this end, we collected ambient noise data while operating temporary seismometers.
2 Data recorded at temporary observatories 2.1 Current status of the KMA observation network Three types of seismometers have been installed at~14 km intervals in the KMA seismic station. These three types, termed surface, borehole A, and borehole B, have been installed at 41, 156, and 85 sites, respectively. Figure 1 shows the distribution of these three seismometer types across the KMA seismic monitoring network. The surface stations are mostly located in areas that are presumed to have a bedrock layer, with borehole stations installed 20 and 100 m below the ground surface. Accelerometers were also installed at 20 m and 100 m below the ground surface. For type-A boreholes, seismometers were only installed 20 m below the ground surface, whereas for type-B boreholes they were installed at both 20 and 100 m below the ground surface. Each type can be easily discerned from the code name of each station, which consists of four   During the period over which the seismic monitoring network has been managed by the KMA, seismometers occasionally malfunctioned. A temporary seismometer (velocimeters or accelerometers) was installed on the ground surface while the malfunctioning sensor was being repaired. Therefore, seismic waves could be recorded at the ground surface and depths of 20 m during the repair period. In this study, those stations at which earthquake data were recorded by temporary seismometers. Table 1 summarizes the selected stations and the earthquake data recorded during the times at which temporary seismometers were installed. These stations are shown in Figure 1. These data were collected based on local magnitudes (M L ) ≥ 2.5 and epicenter distances <150 km. The YOCB station is close to these two locations and recorded the highest number of events during the period of temporary seismometer operation (n = 31). The IMWB and HAMB stations recorded fewer events during this period (n = 10 and n = 16, respectively). When seismic wave monitoring records are available for both borehole and surface stations, the amplification of the soil layer can be easily identified. The amplification can be calculated as the ratio of the peak ground acceleration (PGA) in both positions:

Analysis of event data at the temporary station
where PGA s and PGA wr are the maximum horizontal accelerations recorded at the surface ground (depth = 0 m) and boreholes (depth = 20 m), respectively. The YOCB station had the lowest peak ground acceleration (PGA) ratio between the rock outcrop and subsurface recordings at 3.85. Table 2 summarizes the average PGA ratios of surface-to-within motions for the earthquakes. In terms of average PGA ratios throughout the overall event case, the YOCB site shows a small value of 4.10, as expected. However, this value is still higher than the maximum SAC used in the Republic of Korea (MOLIT, 2018;Cho et al., 2022). The disparity in the recorded data for both indicates that the amplification table for seismic design (which is based on rock outcrop measurements) is unsuitable for estimating horizontal seismic waves at ground surface vibrations using borehole stations. This is based on the theory of wave propagation, as shown in Figure 3. Figure 3 shows that the measurement of the borehole sensor affects both the incidence wave in rock and the transit wave of soil layers over rock. Therefore, the application of the SAC method in borehole stations is barely suitable due to the differences between the borehole and rock outcrop sensors.
If records from both the ground surface and within a rock were available, the site effect in the frequency domain from an earthquake could be calculated. The site effect indicated the amplification ratio of both the ground surface and within a rock. We called it the transfer function (TF), which is calculated as follows: where H(ω) is the TF. Amp s (ω) and Amp wr (ω) are the amplification in the frequency domain at the surface ground (depth = 0 m) and boreholes (depth = 20 m), respectively. Figure 4 illustrates the average TF of the layer, from the depth of the borehole sensor to the ground surface. The three stations used for the analysis had a temporary sensor installed, which was the same type as the one on the ground surface due to the failure of the 100 m depth velocimeter sensor. The surface sensors were velocimeters, whereas those in the boreholes were accelerometers. Therefore, the data recorded by the surface and borehole recordings were converted to m·s -2 , in line with physical quantities. To characterize the signal in the frequency domain, a fast Fourier transform (FFT) was performed, and the smoothing technique of Konno and Ohmachi (Konno and Ohmachi, 1998) was applied. The value of b in the method was set to 100. The TF was calculated as the spectral ratio between surface and borehole signals.
The fundamental frequency (f 0 ) at IMWB and HAMB stations was 7 and 10 Hz, respectively. The f 0 is associated with a velocity contrast in a shallow soil layer (Zhu et al., 2020), and is found first mode of peak in the horizontal-to-vertical spectral ratio (HVSR) (Kwak and Seyhan, 2020). The time-averaged shear-wave velocity to a depth of 30 m (Vs30) was estimated for each of these two stations using the empirical formula suggested by Hassani and Atkinson (2016), and the calculated values were 540 and 676 m s −1 , respectively. These estimates correspond to "very dense soil" and "soft rock," respectively, based on the National Earthquake Hazards Reduction Program Site Classification Standards. In the same way,

FIGURE 3
Schematic of difference between transfer function (TF) and site amplification coefficient (SAC).
Frontiers in Earth Science frontiersin.org  Frontiers in Earth Science frontiersin.org 06 the f 0 at the YOCB station is 15 Hz, with an estimated Vs30 of 872 m s −1 , which his classified as "rock".

Collection of field data
Ambient noise data were obtained by segmenting data collected on a day without an event. Figure 5 shows photographs of the 3 stations and their locations. The area of HAMB comprises farmland, and a stream flows in the vicinity of the seismic station. The IMWB station is located on a mountain and is close to a sewage treatment plant. The YOCB station is inside an Automatic Weather Station, on a plain. The three stations have a 100 Hz sampling rate and a Q330HRS recorder (Kinemetrics). IMWB has an STS-2.5 sensor, while HAMB and YOCB both have STS-2.0 sensors (Kinemetrics).
The 24-h waveform at each station was divided into 48 timewindows of 30 min, as artificial noise levels are dependent on the timing of human activity (Roy et al., 2021). Figure 6 compares the ambient noise data recorded by the borehole and surface sensors. The monitoring time for the compared data was 22:10:00-22:11:40 (KST), when there was generally relatively little human activity. Thus, the influence of artificial noise was presumed to be negligible. The noise level in EW direction at the IMWB station was lower than that of the vertical component. As IMWB station was located in the middle of the mountain, a vertical component occurred due to topographical effects (Del Gaudio et al., 2018). We calculated power spectral density (PSD), which are shown in Supplementary Figure S1.

Data analysis 4.1 Data preprocessing
Data analysis conducted in this study used vibrations recorded by the borehole and surface seismometers at identical time points. Two analyses were performed using the surface-to-borehole vibration ratio. First, the surface-to-borehole ambient noise ratio was calculated. The Fourier amplitude spectrum (FAS) was used to compare the ambient noise; here, only the horizontal components of the noise were considered. The FAS procedure suggested by Ahn et al. (2021) was used as follows: Secondly, to incorporate the influences of the vertical components, HVSR was applied following the FAS. Many studies have shown that the fundamental frequency (f 0 ) of a site can be obtained from the predominant frequency of the ambient vibration HVSR (Sylvette et al., 2006;Guillier et al., 2008;Nagashima et al., 2014;Hassani and Atkinson, 2016;Oubaiche et al., 2016;Kwak and Seyhan, 2018;Ahn et al., 2021). In fact, a Rayleigh wave with an elliptical waveform leads to the disappearance of the vertical amplitude at f 0 , and the Airy phase of a Love wave leads to the collision of horizontal energy at f 0 . Using these two phenomena, the reliability of the f 0 evaluation method has been theoretically and empirically verified for a simple 1D geological structure, assuming that the peak of the H/V ratio approximates f 0 and there is a strong impedance contrast (Sylvette et al., 2006). The interpretation of the HVSR curve is complicated by a lack of understanding regarding the precise composition of the microtremor wavefield, however, the characteristics of the soil layer from the amplitude and frequency could be estimated (Molnar et al., 2022).

Ambient noise ratio of surface to within motion
Ambient noise is a useful variable because it can easily provide site information about the soil layer. We used the HVSR method, which has the advantage of estimating the site effect without the influence of the source and propagation paths on the response spectrum (Xu and Wang, 2021). Based on this principle, the difference in noise at each depth was analyzed at the same point. Both the HVSR and the differences between components were compared. Figure 7 illustrates the HVSR of the borehole and surface sensors. The amplitude ratio of HVSR depends on the difference between the horizontal and vertical components. Thus, the HVSR of within-motion at the borehole seismometer is small relative to the ground surface motion. The maximum amplitude at the HAMB station occurs at approximately 10 Hz, which is similar to the results obtained from the TF based on the event. However, the HVSR at the ground surface of IMWB and YOCB stations has a low amplitude, and the dominant frequency is not clearly shown. The vertical component in the IMWB was greater than that at the other stations ( Figure 6). Ultimately, it is suggested that the HVSR is relatively low because of the high vertical component noise. The HVSR of YOCB is due to a thin soil layer, as previously analyzed.
The left column of Figure 8 shows the surface-to-withinmotion ratio (SWMR) from the borehole to the surface sensors based on ambient noise. The right column of Figure 8 illustrates the SWMR from the HVSR, termed SWMR HVSR . At the HAMB station, the maximum ratio occurred at a frequency of approximately 10 Hz in both analyses. We found that the resonance frequency that is not visible in the HVSR at IMWB and YOCB stations is slightly clearer in SWMR. The resonance frequency at the IMWB station was determined to be within the range 6.5-7.5 Hz from the NS component. At this time, IMWB seems to show confirmed the amplification of a specific period even in the vertical component due to the topographical effect (Massa et al., 2014;Wang and Sun, 2019). The resonance frequency at the YOCB station was determined to be approximately 16 Hz, based on the event, ambient noise, and HVSR data. The amplitudes of both SWMR and HVSR increased with increasing frequency, while SWMR HVSR did not. This response may have resulted from the overlapping effect of the ambient noise propagation reflected by the bedrock, but we did not identify the cause. The HVSR at the YOCB station is influenced by the hard bedrock layer at shallow depth, overlain by a shallow soil layer. This case has similar characteristics to those reported by Ahn et al. (2021). Such characteristics are common when the soil layer is thin and the Frontiers in Earth Science frontiersin.org two layers have a high impedance ratio (Ghofrani and Atkinson, 2014;Kwak and Seyhan, 2020). Overall, the noise-based SWMR (SWMR n ) was higher than or similar to the estimated TF from the event. The gap between the SWMR n and TF is big in the high-frequency domain. This indicates that the ground surface has higher artificial noise levels of unknown high frequencies than the underground. As the high frequencies are steeply attenuated (Peng et al., 2019;Stanko et al., 2020), there is a weak propagation of artificial noise of unknown high frequencies in the underground. For this reason, artificial noise of small amplitude may not be recorded on seismometers 20 m below the ground surface. Additionally, the SWMR HVSR was found to be similar to or lower than the TF identified during the event.

Proposed model
In the previous analysis, we confirmed that the noise-based HVSR and SWMR are similar to the event TF in the frequency domain. However, the horizontal amplitude of the ratio of both the HVSR and SWMR depends on the topographical characteristics of the geology (Del Gaudio et al., 2018;Ahn et al., 2021). Therefore, we need an estimation method for calculating TF using SWMR or HVSR.
We proposed a method to correct the horizontal component ratio of between surface and within-rock by the vertical component ratio. The proposed equation is as follow: where H_Amp and V_Amp are each of horizontal amplitude and vertical amplitude. The subscripts "s" and "wr" indicate the location of the sensor, in either the ground surface or borehole (within rock), respectively. The proposed method combines the noise ratio based on the surface-to-borehole ratio of ambient noise and HVSR results of the surface data. Figure 9 shows noise-based TF methods, compared with the eventbased TF. The proposed model is similar to the TF based on the event and shows excellent performance in the amplitude ratio of the predominant frequency. However, there was a difference in the characteristics at high frequency (20 Hz over). We verified the reliability of the proposed model using recording data in the next chapter.
6 Results and discussion

Comparison of numerical method and proposed model
To verify the estimated surface motion based on TF, these values were compared with the ground surface records. First, the borehole record was calculated with the FFT and multiplied by the TF. Then, the calculated frequency value was converted into acceleration time history using inverse FFT again. This allowed us to obtain the ground surface motion.
In Figure 10 , we show the horizontal components of the surface and synthetic seismograms for each station.
Overall, the estimated waveforms were similar to those of the recorded data at the ground surface. Although there were differences between the estimation approach and both methods for ground surface motion, and minor. More specifically, the PGA obtained using the surface motion at the HAMB station was 0.057 m·s −2 , whereas that obtained from the proposed model was 0.075 m·s −2 (a difference of 0.018 m·s −2 ). For the IMWB station, the PGA based on the surface monitoring data was 0.124 m·s −2 . The PGA estimated using the proposed model was 0.143 m·s −2 . The PGA of the YOCB station at the ground surface was 0.037 m·s −2 , whereas the estimated PGA was 0.024 m·s −2 (a difference of 0.013 m·s −2 ). The PGA gap shows a slight deviation from the recorded surface data. Figure 11 shows the 5% damped acceleration response spectra of the input and surface ground motions (proposed model and recorded) for the case shown in Figure 2. Overall, the spectral accelerations of ground motions estimated by proposed model were similar to the surface recorded data at all stations. More specifically, the proposed IMWB and YOCB models underestimated the spectral accelerations in the natural frequency region compared to the surface record. Although the difference between the two results is minor, we decided to use the proposed model because it is more reliable than SAC.
As our aim to develop a technology that will improve the accuracy of the KMA's seismic intensity information service, we analyzed the accuracy of the calculated MMI. Figure 12 shows the MMI obtained by applying the proposed model and SAC of the Borcherdt model to the borehole data, as well as the MMI obtained using the surface data. The Borcherdt model is currently being applied to the seismic intensity maps at the KMA. Figure 12A shows MMI obtained by applying the Borcherdt model for the borehole record. To obtain it, the amplification coefficient of Table 3 was applied for each period after the borehole sensor record was frequencyconverted. Different MMI conversion equations from the developed KMA model were applied as follows: where PGA was applied and expressed a unit in cm·s −2 . SAC correction data are expressed as MMI SAC and surface data are expressed as MMI Observed . The MMI based on the SAC of the Borcherdt model tended to slightly underestimate that of the surface data. However, the MMI based on the proposed model displayed a 1:1 relationship with the surface data, thereby verifying the outstanding performance of the proposed model.

Verification
We analyzed noise with three stations and proposed the method of a transfer function to predict. However, the results from the three cases may not be able to represent all stations. Therefore, we additionally searched stations that had one or more events during the period of the temporary surface seismometer operation and used them for a verification case. Although we know the estimated value of Vs30 (Ahn et al., 2021;Cho et al., 2022), there are no shear wave profiles available for the 21 stations used in the verification. For these stations, we collected ambient noise records on both the ground surface and borehole within rocks. We calculated the SWMR and HVSR, which are shown in Supplementary Figure S2. The transfer function was calculated based on the proposed method, as shown in Supplementary Figure S3.

FIGURE 8
Comparison of surface-to-within motion ratio (SWMR) based on ambient noise. Panels on the left represent the horizontal spectral ratios (H surface / H borehole ); center panels are the vertical spectral ratios (V surface /V borehole ); and right panels are the HVSR ratio (HVSR surface /HVSR borehole ).  The results of verifying the transfer function based on the proposed method for 45 events in 21 station records are shown in Figure 13. The seismic intensity was calculated by Eqs 4, 5, and the first decimal place was rounded up and expressed in Roman letters. In Figure 13, outside of the square the case of under-or over-predicting the seismic intensity is reported. In the Borcherdt model-based KMA process, the accuracy

Automatic computation of seismic intensity
Seismic stations provide seismic intensity information that represents seismic intensity value at the local site and city. If KMA provided the borehole sensor records within motion without correction, the intensity may be sorely underestimated. Therefore, we recommend that the proposed method be applied to each station.
We have designed an automatic computation for seismic intensity, which follows the procedure outlined below: (1) Extract a record at the 20 m borehole accelerometer sensor after the earthquake (time window default: 300 s) (2) Apply a noise filtering to the waveform using Choi et al. (2019) method (band-pass filter default: 0.1-50 Hz) (3) Apply FFT to the event waveform, and multiply the TF (TF is provided as Supplementary Material) (4) To extract the waveform, apply inverse FFT to the results of 3 step. (5) Calculate PGA of the waveform, calculated as MMI (using Eqs 4, 5).

Conclusion
The aim of this study was to estimate horizontal direction waves recorded by in-borehole seismometers in order to accurately estimate seismic intensity at the ground surface. The proposed model is based on ambient noise data collected from both surface and borehole sensors.
To this end, the fundamental frequency of the soil layer was determined based on the TF values of the event, ambient noise, and HVSR without performing ground surveys. At locations with both borehole and surface sensors, the natural frequency of the soil layer can be calculated via a surface-to-within-motion ratio based on the HVSR.
At high frequencies, the SWMR was generally high with respect to ambient noise and relatively low for SWMR HVSR . However, the noise-based TF calculated using the proposed model exhibited a trend similar to the TF of the event analysis. The TF values obtained via event analysis mostly came from intact ground data, which allowed for the analysis of linear behavior. Thus, the model proposed in this study did not consider non-linear amplification; further studies should therefore address this consideration.
The proposed approach requires a temporary seismometer to be installed on the surface for 1 day. The TF can be calculated using the obtained ambient noise. If we perform frequency domain response analysis based on the TF, we can estimate the acceleration time history on the surface. This approach has a significant advantage of not requiring additional geotechnical surveys. However, the present study only analyzed three stations, and so the approach could be further improved through additional research.
The validation of the proposed model showed that it can be applied to borehole data to estimate ground surface vibrations. Although the exact waveform generated by the model did not always exactly recreate observational data, the model delivered PGA values with lower errors than those obtained using the site coefficient. Therefore, this model represents a significant improvement in the accuracy of seismic intensity estimation (for which the site coefficient is currently employed). Using the TF developed in this study could enable reliable seismic intensity estimations to be obtained while operating borehole data; this would be particularly advantageous at stations that lack surface seismometers such as KMA.
Finally, the automatic calculation process based on the proposed method can be applied in real-time. Most earthquake information is provided in a brief time immediately after an event. At this time, the accuracy of the information affects the reliability of the institution or the government and future research. Therefore, the provided automated real-time MMI processing will be very useful.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.