We collected 2022 vertical-component digital seismograms from 155 seismic events recorded at 93 stations in and around the Korean Peninsula. The seismic events included 146 natural earthquakes between May 2010 and May 2022, six North Korean nuclear tests, and three chemical explosions that occurred in August 1998 (Zhang et al., 2002; Yang et al., 2021). These seismic events occurred within the crust at focal depths shallower than the Moho discontinuity from CRUST1.0 (Laske et al., 2012). Natural earthquakes with moderate magnitudes between m _b 3.5 and 5.3 were selected to ensure high SNRs and to avoid complex rupture effects. The epicentral distances are between 150 and 1700 km to observe the Lg-wave phase clearly (e.g., Zhao et al., 2010; 2013). We visually inspected individual traces to remove low-quality data that were saturated or noisy or had incorrect timing, possibly due to off-point records, low magnitudes, or the superposition of multiple events. Although the selection process reduces the amount of available data, a reasonable dataset is obtained following this process. The locations of the stations and events used in this study are shown in Figure 1, and their parameters are listed in Supplementary Tables S1, S2.

Following Zhao et al. (2013), we processed the seismic data to extract Lg-wave amplitude spectra. After removing the trends, means, and instrument responses from the raw vertical-component seismograms, we scanned the most energetic waveforms within a group velocity window of 0.6 km/s between 3.7 and 2.8 km/s for Lg amplitude measurements (Figure 2) (Zhao and Xie, 2016; Zhang et al., 2022). Subsequently, we captured the time series of pre-P noise and pre-Lg noise with a time length window equal to the Lg waveform and calculated the waveform energy. The fast Fourier transform was used to calculate the amplitude spectrum of the Lg waves and noises at 66 discrete frequencies log-evenly distributed between 0.05 and 20.0 Hz. Pre-P noise and pre-Lg phase noise were used for both data quality control and correction of the Lg-wave energy (e.g., Luo et al., 2021). Lg waves with an SNR to pre-P noise of less than 2.0 were removed to ensure Lg data quality. Furthermore, we also set a pre-Lg SNR threshold of 1.0 to remove data that were possibly dominated by Sn coda (Zhang et al., 2022). Then, the amplitude spectrum of the Lg wave can be obtained by using A S 2 = A O 2 − A N 2 , where A S , A O , and A N are the amplitude spectra at different frequencies of the true signal, the observed data, and noise, respectively. After data screening and denoising, the Lg-wave spectral dataset can be obtained for Lg-wave attenuation tomography.

The Lg-wave amplitude spectrum recorded by station i from event k at frequency f can be expressed as follows (Xie and Mitchell, 1990) A k i f , ∆ = S k f G k i ∆ Γ k i ∆ , f P i f r k i f , (1) where ∆ is the epicentral distance, S k f is the source term, G k i ∆ = ∆ 0 ∆ k i − 1 / 2 is the geometrical spreading function with a reference distance of ∆ 0 = 100   k m (Street et al., 1975), Γ k i ∆ , f is the attenuation term, P i f is the site response and r k i f is the random effects of minor factors along the propagation path and computational errors.

The attenuation factor can be expressed as follows (Zhao et al., 2013) Γ k i Δ , f = ⁡ exp − π f v ∫ k i d s Q x , y , f , (2) where v represents the Lg-wave group velocity, ∫ k i d s is the integral along the great circle from event k to station i , and Q x , y , f is the Lg-wave quality factor related to the surface location coordinates x , y .

When two stations i and j record the same event k and the locations of the stations and event are approximately aligned on a great circle, the two-station Lg-wave amplitude ratio from stations i and j can be calculated as follows (Zhao et al., 2013) A i j = A k j A k i ≈ ∆ k j ∆ k i − 1 2 ∙ ⁡ exp − π f v ∫ i j d s Q x , y , f ∙ P j P i ∙ r k j r k i , (3) where A k i and A k j are the observed amplitudes at stations i and j from a single event k , ∆ k i and ∆ k j are the epicentral distances from k to i and j , respectively, ∫ i j d s is the integral along the great circle from i to j , P i and P j are the site responses at i and j , respectively, and r k i and r k j are random errors along the ray paths from k to i and j , respectively. Compared with the single-station data in Eq. 1, the two-station data shown in Eq. 3 effectively eliminate the compromise error between the attenuation and source.

To establish an inversion system for Lg-wave Q tomography, we apply the natural logarithm to Eqs 1, 2 based on perturbation theory (e.g., Zhao et al., 2010; 2013). By neglecting the random effects along the propagation path, we assume that r f = 1 ; then, we have ln A k i f , ∆ = ⁡ ln S k f + ⁡ ln S i ∆ − π f v ∫ k i d s Q x , y , f + ⁡ ln   P i f , (4)

The terms Q x , y , f , ln S f and ln   P f can be separated into a background value and a perturbation value (Zhao et al., 2013): 1 Q x , y , f ≈ 1 Q 0 x , y , f − δ Q x , y , f Q 0 x , y , f 2 , (5) ln S k f = ⁡ ln S k 0 f + δ ⁡ ln S k f , (6) ln P i f = ⁡ ln P i 0 f + δ ⁡ ln P i f , (7) where Q 0 , S 0 and P 0 are the background values of the Lg-wave quality factor, source spectrum and site response spectrum, respectively, for beginning the inversion. By substituting Eqs 5-7 into Eq. 4, we have ln A k i f , ∆   = ⁡ ln A k i 0 f , ∆ + δ ⁡ ln S k f + π f v ∫ k i δ Q x , y , f Q 0 x , y , f 2 d s + δ ⁡ ln P i f , (8) where ln A k i 0 f , ∆ = ⁡ ln S k 0 f + ⁡ ln G k i ∆ − π f v ∫ k i d s Q 0 x , y , f + ⁡ ln P i 0 f . (9)

Then, we obtain the amplitude spectrum residual by taking Eq. 8 minus Eq. 9 as follows: δ ⁡ ln A k i f , ∆ = δ ⁡ ln S k f + π f v ∫ k i δ Q x , y , f Q 0 x , y , f 2 d s + δ ⁡ ln P i f . (10)

Therefore, the perturbations in the attenuation, source, and site response are related to the amplitude residual. The amplitude residual δ ⁡ ln A k i f , ∆ for event k recorded by station i at frequency f is denoted as h ∼ k i f . It can be distributed to the path based on the following mesh discretization: h ∼ k i f = ∑ n = 1 N a i n · δ Q n + e k · δ ⁡ ln S k f + u i · δ ⁡ ln P i f , (11) where n is the index of a grid point, N is the total number of grids of a ray path, a i n = − π f v D n Q 0 x n , y n , f 2 , x n , y n and D n are the coordinates and the length of the ray path in grid n , and e k and u i are coefficients for event k and station i , respectively, with e k = u i = 1 for single-station data. Then, we have the linear matrix equation of the Lg-wave Q perturbation for single-station data H s = A s · δ Q + E s · δ S + U s · δ P , (12) where H s is a vector composed of Lg amplitude spectra residuals, δ Q is a vector composed of the perturbations of the Q models, matrix A s is composed of elements a i n and sets up the relationship between Q perturbations and the observed Lg-wave spectra, δ S is a vector composed of the perturbations of source terms, matrix E s sets up the relationship between the source perturbation and the observed Lg-wave spectra (Zhao et al., 2010), δ P is a vector composed of the perturbations of site response terms, and matrix U s sets up the relationship between the site response perturbations and the observed Lg-wave spectra (Zhao and Mousavi, 2018).

For two-station data, since the source term is eliminated by taking the spectral ratios in Eq. 3, the similar linear matrix equation of the Lg-wave Q perturbation is H d = A d · δ Q + U d · δ P . (13)

In general, it is assumed that ∑ n = 1 N δ ⁡ ln P i f = 0 when the stations are evenly distributed (Ottemöller et al., 2002; Zhao and Xie, 2016; Zhao and Mousavi, 2018); thus, the site response terms δ P are commonly ignored in the inversion of single- and two-station data. By combining Eqs 12, 13, we have H s H d = A s A d · δ Q + E s 0 · δ S . (14)

Eq. 14 forms a joint inversion problem for perturbation vectors δ Q and δ S . The regional average Q obtained from the two-station data is used as the initial model. The inversion can be solved by the least-squares QR method, which includes regularization, damping and smoothing (Paige and Saunders, 1982). The current Lg-wave Q correction is used in the next iteration until satisfactory convergence is obtained, and there are 250 iterations at each frequency from 0.05 to 20 Hz. Following inversion, the amplitude residuals are closer to the Gaussian distribution, and the root mean square values of the total residuals at all 66 frequencies are significantly reduced (Figure 3).

The checkerboard method was used for the resolution analysis (Zelt, 1998; Morgan et al., 2002; Zhao et al., 2013). To create a checkerboard-shaped Q_Lg model, ±7% logarithmic perturbations are superimposed on a constant background Q_Lg (Zhao et al., 2013). We generated synthetic Lg spectral amplitudes according to actual ray paths (Figure 4A) and added 10% random noise to form a tomographic dataset (e.g., Zhang et al., 2022), where the source functions adopted inverted source spectra. Both single- and two-station synthetic data were jointly used to reconstruct the Q_Lg model (e.g., Zhang et al., 2022). Figure 4 shows the ray path coverage and reconstructed Q perturbation model at 1.0 Hz. Note that the resolution is significantly improved compared with that of previous studies (Zhao et al., 2010; 2013). The workflow for data collection, pre-processing, tomography, and verifications can be referred to Zhao et al. (2022).

Based on the inversion procedure described in the previous section, we obtained the Q_Lg model at 66 individual frequencies between 0.05 and 20 Hz. Figure 5 shows the Q_Lg map coverages at 0.5, 1.0 and 2.0 Hz. The lateral Q_Lg variations are consistent with the regional tectonic conditions. In the Korean Peninsula region, the Q_Lg distribution is characterized by high values for the Nangrim Massif (NM) in the north and the Yeongnam Massif (YM) in the south and low values for the Gyeonggi Massif (GM) in the middle, especially at a frequency of 1.0 Hz (Figure 5B). Our results are consistent with those of previous studies but at relatively higher resolutions, and they cover the entire Korean Peninsula (Zhao et al., 2010; 2013; Pasyanos et al., 2012; Ranasinghe et al., 2014).

We investigated the regional variations and frequency dependence of Q_Lg in three massifs to characterize the Lg-wave attenuation for different geological formations on the Korean Peninsula. For example, the scattered distribution of observed Q_Lg values is shown for the Gyeonggi Massif (Figure 6A), whereas the average Q_Lg versus frequency is summarized for all three massifs in the Korean Peninsula (Figure 6B). The average Q_Lg values increase with increasing frequency for an individual massif below 0.3 Hz and above 3.0 Hz. However, the Q_Lg variation in the Gyeonggi Massif (GM) differs from that in the other two massifs between 0.3 and 3.0 Hz.

Figure 7 shows a comparison of the resulting broadband Q_Lg with the surface topography, crustal thickness, and seismicity. The Lg attenuation variations indicate changes in the crustal waveguide and crustal physical properties across the three massifs. The boundary between the NM and GM corresponds to low Q_Lg values of approximately 1.0 Hz (Figures 7C, D). Q_Lg usually decreases with decreasing crustal thickness (Zhao et al., 2010; 2013). Numerical simulation studies have shown that the seismic Lg wave attenuates with crustal thinning along the waveguide (e.g., Hong et al., 2008). However, the abnormally low Q_Lg for the GM does not correspond to a significant Moho depth increase in the Korean Peninsula (Figure 7B). Therefore, the Lg attenuation is not strongly affected by a smooth anomaly in the Moho depth (Campillo, 1987). The thickness variation in the crustal waveguide is not the main factor affecting the Q_Lg on the Korean Peninsula. Strong seismicity can be observed in the uppermost mantle beneath the low-Q_Lg crust at approximately 38.5°N (Figure 7B). Therefore, the low-Q_Lg region, located at the boundary between the NM and the GM, can be attributed to complex tectonic sutrue between two ancient plates (the North China Craton and the South China Block), where the crustal structure has been strongly influenced by tectonic processes, including the extension of the Sulu collision Belt (Chough et al., 2000; Zhai et al., 2019).

During the joint inversion, the Lg-wave source excitation spectral amplitudes at 66 discrete frequencies are calculated for all events. The source parameters are obtained by fitting the resulting Lg excitation spectra (Zhao et al., 2010; 2013; He et al., 2020). We calculated the scalar seismic moment M 0 , the corner frequency f c , and the high-frequency falloff rate n by fitting the Lg-wave excitation spectrum with the ω − n source model (Brune, 1970; Street et al., 1975; Sereno Jr et al., 1988). The source term in Eq 1 is expressed as follows: S f = M 0 4 π ρ v s 3 1 + f f c n , (15) where ρ and v s are the average density and shear-wave velocity in the crust, respectively, with values of 2.7 g/cm³ and 3.5 km/s for Northeast China (Zhao et al., 2010). Figure 8 shows the best-fit source models with solid colored lines, and the shaded areas represent their standard deviations. Then, M 0 , f c , and n can be determined by minimizing the L2 norm of the residuals between the theoretical source function and the network-determined source spectral data (Figure 8). The M 0 for NKT1-6 increased successively. According to the Mueller–Murphy model (Mueller and Murphy, 1971), the source corner frequency ( f c ) is predicted to decrease with increasing yield and increase with increasing source depth. According to the source model parameters, the explosion yields increase for NKT1-6; however, the corner frequencies, represented by f c , are not consistently reduced due to differences in burial depth. Figure 8 shows the source excitation spectra inverted from the observed Lg-wave spectra for NKT1-6. Figure 9 provides a comparison between the retrieved spectra and the synthetic spectra. Except for NKT3, the source spectra gradually increased from NKT1 to NKT6 at all frequencies. The NKT3 source spectra are larger than those of NKT5 at high frequencies (above 0.5 Hz); thus, NKT3 is larger according to the previously estimated m _b (Lg) based on an Lg-wave amplitude of approximately 1.0 Hz (Zhao et al., 2014), which is inconsistent with the order of m _b (P) from the USGS National Earthquake Information Center (NEIC) and the Comprehensive Test Ban Treaty Organization (CTBTO).

Figure 10 shows a comparison between the seismic moment M 0 for NKT1-6 and the results from other studies. Chiang et al. (2018) used a time-domain waveform inversion to calculate the full moment tensor from regional stations in China, South Korea, and Japan with the MDJ2 1D-layered Earth model (Ford et al., 2009) to calculate Green’s functions. Alvizuri and Tape (2018) applied a grid search and minimized the misfit function between observed and synthetic waveforms to determine the full moment tensor, and they used the MDJ2 1D-layered Earth model for calculations. The period band of analysis by Chiang et al. (2018) was typically 20–50 s, with shorter periods used for NKT1. The M 0 results in this study are consistent with those of Chiang et al. (2018) for NKT2-6, with a large deviation for NKT1. In the low-frequency (<0.05 Hz) band, crustal attenuation variations could be ignored. Therefore, the M 0 consistency in NKT2-6 verifies that the Lg source spectra accurately remove the path attenuation effect in this study, and inconsistency for NKT1 might be related to poor SNR at low frequency.

Pasyanos (2022) estimated the yields of NKT1-6 by utilizing the seismic moment function of Denny and Johnson (1991), in which the moment is proportional to the yield. Replacing the units of yield from kilotons to joules, the moment-to-yield ratio ( M 0 / W ) in units of N⸱m/J is obtained as follows: M 0 / W = 3.76 × 10 − 3 V P 2 V S − 1.1544 ρ 0.5615 z − 0.4385 10 − 0.0344 G P , (16) where the material properties ( V P , V S , ρ , and G P ) indicate the P-wave velocity, S-wave velocity, density, and gas porosity, respectively. The material properties ( V P , V S , ρ and G P ) used in this study are 5500 m/s, 3175 m/s, 2550 kg/m³, and 0.02%, respectively, for granite at the NKT site (Stevens and Day, 1985). By using differential travel times from Pn and Pg waves, the burial depths (z) of NKT1-6 were determined to be 330, 540, 506, 468, 521, and 570 m (Yang et al., 2021). With the emplacement conditions and burial depths provided above, the yields can be estimated by dividing the seismic moment by the right side of Eq 16, and the values are 4.6, 8.5, 19.9, 20.9, 24.7, and 337.4 kt (Figure 11).

Figure 11 shows a comparison between the yield estimation results for NKT1-6 above and several other estimates based on (1) regional waveform envelopes (Pasyanos and Myers, 2018), (2) the seismic moment using the formula of Pasyanos (2022) with M 0 from Chiang et al. (2018), (3) the intercorrelation procedure, which applies waveform equalization to teleseismic P and regional Pn seismograms (Voytan et al., 2019), (4) high-frequency (>4 Hz) filtered P waves (Voytan et al., 2019), (5) source spectral ratios of narrow-band regional body-wave coda waveform envelopes (Delbridge et al., 2023), (6) NEIC teleseismic m _b using the empirical relationship of Bowers et al. (2001), followed by a depth correction using the equation of Patton and Taylor (2011) with the burial depths determined by Pn and Pg differential travel times (Yang et al., 2021), and (7) regional m _b (Lg) (Xie and Zhao, 2018). The yield estimations in this study are highly consistent with teleseismic m _b -derived yields. The body wave magnitude (m _b ) based on teleseismic phases is rarely influenced by the crustal structure along the ray path. The purpose of using the Lg-wave to calculate the seismic moment in this study is to determine whether Q_Lg model can be used to eliminate the attenuation effect. The results show that after effectively removing the attenuation effect, the moment rather than the magnitude of the regional seismic Lg agreed with the teleseismic m _b yield estimation.