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 letters, with the first three indicating the name of the region and the last indicating the sensor type (A or B).

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).

Figure 2 shows the accelerograms for the horizontal component borehole and surface at the three aforementioned stations. Figures 2A, B displays the seismic record from the HAMB station, showing an M _L 5.8 earthquake that occurred at 11:32:54 on 09/12/2016 (UTC), 131 km from the station. Figures 2C–F show seismic recordings from IMWB and YOCB stations, respectively, for an M _L 5.4 earthquake that occurred at 05:29:31 on 15/11/2017, at distances of 125 and 40 km, respectively, from the earthquake epicenter.

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: A m p l i f i c a t i o n = P G A   r a t i o = P G A s P G A w r (1) 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: H ω = A m p s ω / A m p w r ω (2) 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 ₀ ) at IMWB and HAMB stations was 7 and 10 Hz, respectively. The f ₀ 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⁻¹, 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, the f ₀ at the YOCB station is 15 Hz, with an estimated Vs30 of 872 m s⁻¹, which his classified as “rock”.

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:

(1) FFT was applied to the horizontal components of the motions.

Secondly, to incorporate the influences of the vertical components, HVSR was applied following the FAS. Many studies have shown that the fundamental frequency (f ₀ ) 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 ₀ , and the Airy phase of a Love wave leads to the collision of horizontal energy at f ₀ . Using these two phenomena, the reliability of the f ₀ 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 ₀ 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 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-within-motion 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 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.

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⁻², whereas that obtained from the proposed model was 0.075 m·s⁻² (a difference of 0.018 m·s⁻²). For the IMWB station, the PGA based on the surface monitoring data was 0.124 m·s⁻². The PGA estimated using the proposed model was 0.143 m·s⁻². The PGA of the YOCB station at the ground surface was 0.037 m·s⁻², whereas the estimated PGA was 0.024 m·s⁻² (a difference of 0.013 m·s⁻²). 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 frequency-converted. Different MMI conversion equations from the developed KMA model were applied as follows: M M I = 1.8976 ∙ log P G A + 1.8365   I ≤ M M I ≤ V (4) M M I = 2.8828 ∙ log P G A + 0.3945   V < M M I ≤ VIII (5) where PGA was applied and expressed a unit in cm·s⁻².

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.

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.

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 of the MMI prediction was 7.3 %, while the accuracy of the MMI calculated by the proposed method was 69.1 %. We confirmed that the accuracy of predicting improved by the proposed model.

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)