Direct Estimation of Explosion Pn Green ’ s Functions Applied to Teleseismic P -wave Intercorrelations for North Korean Nuclear Tests

Seismic yields and explosion source time functions for the six declared underground nuclear tests at the North Korean test site are estimated using teleseismic waveform equalization (intercorrelation) incorporating improved source structure inferred from regional distance Pn phases. Explosion source time functions estimated by previous modeling efforts are deconvolved from the broadband Pn observations to extract direct estimates of the Pn Green ’ s functions. Very consistent signals are found for station MDJ (to the north) for all but the 2006 event, and a 1D layered elastic model is determined that matches the corresponding pPn signals for the larger events. The much simpler Pn Green ’ s functions found for INCN (to the south) are well modeled using a half-space structure. The 1D structures are incorporated into intercorrelations for teleseismic P observations. Very similar results are found using the layered MDJ model for all azimuths, or when azimuthally varying the source velocity model to have a half-space structure for stations to the south, and/or using a half-space structure for the 2006 event at all azimuths. The relative yield estimates from intercorrelation are thus stable with respect to specific 1D source structures; the full effects of 3D elastic and nonelastic structures are yet to be considered.


Introduction
Reliable estimation of nuclear explosion yields for test sites for which we do not have data from a direct calibration experiment presents the difficult problem of calibration of seismological properties of the crust and upper mantle around the site. Soviet test sites were calibrated by detailed waveform analysis of teleseismic recordings of underground tests, resulting in m b -yield relationships that were subsequently validated by the Joint Verification Experiment direct calibration. A similar challenge now exists for empirical calibration of the six nuclear explosions conducted by the Democratic People's Republic of Korea (DPRK: North Korea). Current methods for yield estimation for this test site largely rely on transporting various magnitude-yield relationships from calibrated test sites to this uncalibrated site. This is very uncertain given the lack of seismicity in North Korea from which to determine seismological properties that affect seismic-wave excitation and anelastic attenuation of the explosion signals.
Regional and teleseismic seismic-wave recordings for the six underground DPRK nuclear tests have now been modeled in many studies (see reviews in Walter andWen, 2018 andVoytan et al., 2019). Although some procedures attempt to determine source structure and/or path effects, most studies make unverified assumptions or use transported relationships that do not account for possible test site bias. Most procedures do provide fairly consistent relative yield estimates from the diverse seismic measures, but absolute yields for a given event still span a large range, even exceeding the classic "within a factor of two" uncertainty long attributed to seismically determined yield estimates. Absolute yield estimation requires a parameterization of the explosion source function, typically involving dependence on emplacement medium, depth of burial, water saturation, and explosive yield (e.g., Mueller and Murphy, 1971;Denny and Johnson, 1991), along with seismic phase path-specific attenuation. It may also be affected by nonlinear damage and spall effects around and above the source, along with any release of tectonic stress, scattering interactions with test site surface topography, and/or path-specific structure. Constraining these contributing factors from remote observations alone is daunting but necessary, given the lack of immediate prospect of direct calibration at the North Korean test site using a known device.
The approach taken in this article is to seek direct empirical constraint on the near-source velocity structure above the nuclear explosions and to use that improved information (compared to the earlier assumption of a half-space) in a relative waveform equalization procedure called intercorrelation (Lay et al., 1984;Lay, 1985;Voytan et al., 2019) that accounts for source depth, near-source structural effects on depth phases, and parameterized explosion source function dependence on depth and yield, while eliminating uncertain remote path effects. Tied to an absolute yield level based on prior estimation of the test site attenuation factor (average t* for teleseismic P waves), yields are estimated for the six nuclear tests at the North Korean test site.

Starting Estimates of Test Site Attenuation Factor and Explosion Source Functions
The test site attenuation factor was obtained by Chaves et al. (2018) by modeling broadband P-wave ground motions for the 3 September 2017 (NK6) nuclear test, assuming a Mueller-Murphy explosion source time function for a granite source medium and Green's functions that parametrically allowed for pP delays times due to inaccurate above-source velocity structure. The trade-offs between explosion yield, burial depth, and pP behavior, and site attenuation factor (t*) were explored and partially resolved by optimizing the simultaneous fit to farfield broadband P-wave amplitudes and waveforms, along with modeling the absolute amplitude of 4 Hz, high-pass filtered P-wave seismograms for the first 0.25 s of ground motion, as this interval precedes the expected arrival of pP. This combined forward modeling analysis provided an absolute yield estimate (230 ± 50 kt) for the 2017 event associated with a preferred average test site attenuation factor of t* = 0.78 ± 0.03 s.
The 4 Hz amplitudes for the five smaller North Korea explosions were then modeled by Voytan et al. (2019) using the average test site t*-value and Mueller-Murphy models, providing preliminary yield estimates and explosion source functions for all of the events. Again, modeling the 4 Hz signals was designed to minimize effects of pP arrivals in the relative size estimates. The 4 Hz was chosen, because it allows isolation of the first cycle of the filtered data from the later pP arrival. Higher frequency high-pass filtered signals are more attenuated and too noisy to measure reliably. The intercorrelation method was next applied to large global data sets of shortperiod P-wave recordings to determine relative burial depths and adjusted yields of the six DPRK tests, tied to the absolute yield modeling for the 2017 event. The analyses in Chaves et al. (2018) and Voytan et al. (2019) assumed simple half-space Green's functions, given that details of the source crustal velocity structure were unknown, with effects of incorrect shallow seismic velocity structure and/or nonelastic behavior in the source region being explored by allowing for variable pP delay times and/or amplitudes in the intercorrelations. The waveform equalization procedure provides fairly robust estimates of relative source strength with explicit correction for yield scaling and burial depth variations, explicitly accounting for common remote path and receiver effects. If source structure effects are stable to a given station between events, those effects on the differential waveforms are likely accounted for as well. The preferred absolute yield estimates from the intercorrelation procedure given by Voytan et al. (2019) (Table 1) are dependent on the estimated test site average t*-value and on the source velocity structure used for wave excitation. The uncertainty ranges are for a 10% range in the waveform misfit measures, as described in Voytan et al. (2019), and do not include an estimate of model uncertainty. Spectral ratios of regional Pn recordings were compared to predictions for ratios of the absolute source models to validate the differential yield estimates.
The analyses in these prior publications provide estimates of the yields and t* bias for the North Korean test site, enabling calibration of the DPRK magnitude-yield relationship (for any self-consistent set of measured magnitude values) for the test site directly, rather than relying on highly problematic transport of magnitude-yield relationships from other test sites. There are still three large uncertainties associated with the waveform modeling. These are (1) the source region velocity structure is not well constrained and can influence both the coupling of the source and the detailed characteristics of depth https://www.seismosoc.org/publications/the-seismic-record/ • DOI: 10.1785/0320210045 The Seismic Record phases; (2) the mountainous topography at the source region may cause enhancement and azimuthal variations in free surface pP reflections for elastic structures, affecting m b and M s magnitudes (e.g., Avants, 2005;Stevens and O'Brien, 2018); and (3) the shallow depth and size of the 2017 event result in extensive damage around the source and nonlinear interaction with the free surface, which probably affects the pP amplitudes and timing for that event (Chaves et al., 2018;Stevens and O'Brien, 2018). Although the modeling and intercorrelation procedures parameterize for first-order complex effects of pP behavior, it is not clear how much bias may persist in the results due to assuming uniform half-space structure. We address that issue by determining a more reliable layered source structure and using it in the intercorrelation procedure.

Source Function Deconvolutions
The source region velocity structure issue in the list of uncertainties given earlier is addressed by directly determining body-wave Green's functions for the six explosions at regional distance stations that have sufficient bandwidth to provide stable estimates. The general consistency of the foregoing deconvolutions motivates an effort to extract the actual Green's function for each signal (rather than canceling it out, as in the above ratios between events). We use the preferred estimates of the yields and burial depths for the six events from Voytan et al. (2019) listed in Table 1 to compute depth-and yield-dependent granite source RVPs for each event. These are deconvolved from the Pn signals for each available recording at MDJ, INCN, and BJT (Fig. 2). The RVP signals are convolved with a constant Q t* attenuation operator of 0.05 or 0.01 s to suppress high-frequency noise in the deconvolutions and to partially account for path attenuation.
For station MDJ, stable signals are found in the first second of the deconvolutions, with a consistent dovetail feature (two negative peaks) in the downswings for the five larger events suggesting similar upgoing depth-phase complexity, whereas station NK1, located in a different mountain, has a simple  Figure 2 indicate stability of the Green's functions at a given station for the five larger events but with strong azimuthal variation. The empirically derived Pn Green's functions can be compared with the simple halfspace Green's functions used for the regional stations in the previous intercorrelation analysis of Voytan et al. (2019), as shown in Figure 3. These do not account for the dovetail seen at MDJ and underestimate the downswing amplitude. In contrast, the predictions for INCN for the purely elastic half-space model (and for a model with slightly delayed pPn arrival times) are in good agreement with the timing of the weak first downswing in the Green's functions. The half-space model thus appears to be adequate for the signals to the south from the source region but not for signals to the north.
An elastic structure that better matches the pPn impulse responses at MDJ was determined by forward modeling of the estimated Green's functions, emphasizing fitting of the dovetail in the downswing of the pPn energy. The goal is to have a model that can be applied to other northerly stations. Figure 4 shows Green's function predictions for a plane-layered model that provides a much-improved fit to the dovetail feature in the waveforms for the five larger events. The forward-modeled predictions of the observed Pn waveforms (using the RVPs for each event) fit the observed data quite well over the first cycle and a half of the signals. Later arrivals are not predicted, as the path clearly has much more complexity than can be represented by a simple 1D source structure but may be well modeled by calibrated path models like those of Dreger et al. (2021).
This local source structure involves a thin, 10 m thick, very low velocity (strongly weathered) layer with P velocity V P 1:5 km=s and S velocity V S 0:866 km=s, overlying The Seismic Record a 250 m thick layer with V P 3:2 km=s and V S 1:848 km=s, superimposed on a diorite and/or granitic source velocity layer with V P 5:5 km=s, V S 3:175 km=s. The shallow lowvelocity layers may correspond to the weathered basalt cap and stratified volcanic deposits on the north side of Mt. Mantap (e.g., Coblentz and Pabian, 2015). The data are very sparse and do not justify comprehensive modeling, so there is substantial uncertainty and nonuniqueness in the velocity model, with at least 10% uncertainty in the velocities for the three horizontal layers in the structure. Given that the real Earth has rough topography and dipping layers in the mountain, we view this model as a first-order structure that improves on a half-space model, but it must be used with caution. A low half-space V P of 2.5 km/s was estimated from deformation above the NK6 source by Chi-Durán et al. (2021), perhaps representing additional effects of damage for NK6. The dovetail feature shared in the MDJ deconvolutions for NK2-NK6, as modeled by a 1D velocity model, arises from the underside reflections from the low-velocity layer (first negative peak) and the free surface (second negative peak). To the extent that 1D modeling is valid, we can be confident that the sources are deeper than the volcanic layer, given that underside reflections from that layer are observed. Three-dimensional modeling may account for the dovetail feature by a scattering mechanism, so overinterpretation of the 1D structure should be avoided. NK1, located in a separate mountain, which is not believed to have a basalt cap, is not well fit by this structure. Also, the source depth used for NK3 appears to be about 100 m too shallow. The goal here is not to attempt to uniquely characterize the near-source Green's function from the handful of signals but to have simplified 1D elastic structures  Figure 3 are computed using a half-space with the granitic source velocities, so it appears that the southern Pn paths do not sample the volcanic layers. Teleseismic pP signals will have steeper incidence angles, so that the most distant signals at all azimuths could sample some of the shallow layering, and for NK6 these may sense damage around the source volume, but closer teleseismic signals with larger take-off angles should sense the azimuthal variation in the shallow elastic velocity structure for all the events.

Intercorrelation with Estimated Green's Functions
The intercorrelation algorithm was modified to utilize multiple 1D velocity structures, selecting different source structures for stations at different azimuths and/or for different events at the same station. Searches were then performed to find the bestexplosion yields for source depths specified by Voytan et al. (2019) based on location in the mountain and assuming 4% tunnel grades. The intercorrelated waveforms for NK4 and NK5, found using the layered source velocity model inferred from MDJ for all stations, are shown in Figure 5. The waveform  equalization, optimized by minimizing the differential waveform power as in Voytan et al. (2019), works very well for the layered crustal model, giving slightly improved fits relative to using half-space Green's functions, with almost identical yield estimate for NK4 (11.3 ± 0.6 kt for the layered structure versus 11.2 ± 0.6 kt for the half-space with minor pP delay). Holding the NK5 yield at 18.8 kt ties the estimates to the teleseismic t* model. Figure 5 indicates the large number of traces used to measure the relative yields. For these events, an elastic source velocity structure is viable, because the damage zone is expected to be small and to have negligible effects on pP. The events are nearby and likely share any azimuthal pattern in Green's functions, so the waveform differential sensitivity is stable relative to source structure.
Intercorrelations for the other events relative to NK5 using the MDJ-based 1D source structure for all the stations also slightly outperform the prior half-space results with delayed pP arrivals. This stability with respect to significant differences in the Green's functions reflects a key strength of the waveform  Figure 5. Equalized P and Pn waveforms for NK4 (Yield 11.3 kt, depth 710 m) and NK5 (Yield 18.8 kt, depth 710 m). In each station pair, the lower trace is the NK4 signal convolved with the source function (RVP convolved with Green's function for the layered source model for MDJ) for NK5, whereas the upper trace is the NK5 signal convolved with the corresponding source function for NK4. The yield for NK5 was held fixed and the value for NK4 was determined by minimizing the amplitude weighted logarithmic average amplitude error defined by Voytan et al. (2019). Peak amplitudes are normalized to unity for all the traces.
https://www.seismosoc.org/publications/the-seismic-record/ • DOI: 10.1785/0320210045 The Seismic Record equalization approach; as long as the observations at a given station share a common Green's function, the details of the Green's function are not critical. The yield estimates are listed in Table 1 (Table 1). One additional case used the same azimuthal selection of source structures for NK5, but the half-space model was used for all stations for NK1, motivated by the distinct character of the deconvolutions at NK1 in Figures 1-3. The yield estimate for NK1 is then 1.8 ± 0.4 kt. The estimates of yields are not strongly affected by the details of the source function differences, even between sources, particularly when the same source structure is used for all the events observed at a given station. These results are all relative to a fixed yield estimate of 18.8 kt for NK5, the specified source depths, and the Mueller-Murphy granite model.
Because intercorrelation equalizes signals from two events at a given station, for which any shared azimuthal pattern in the explosion Green's functions for similar source depths is common to both the signals, it is likely that extensive computation of Green's functions for all possible source locations in the 3D structure to a large number of teleseismic and regional stations will actually not change the results very much. Simple tests of imposed azimuthally varying pP time for parameterized half-space Green's functions prove very stable as long as the sources share the same azimuthal patterns. Nonetheless, it is clearly desirable to extend the intercorrelation framework to 3D Green's functions in the future. Doing so usefully will require information on the subsurface structure, as sampled here using Pn deconvolutions, rather than simply assuming topography on a half-space. Accounting for nonelastic effects is even more challenging, as this is very yield-dependent, and detailed material properties around the sources are required (and not available), but it appears that only NK6 is likely significantly affected by such effects on pP.

Conclusions
Waveform modeling and waveform equalization procedures provide estimates of explosion yields for the six DPRK test site events based on seismic-wave observations combined with an explosion source model parameterization for a hard-rock site. Interevent deconvolutions of Pn phases for nearby stations demonstrate stable relative behavior of the signals at a given station, with the 2006 test having some differences, and there are strong variations with azimuth for stations MDJ and INCN. Estimates of the explosion impulse responses are made by deconvolution of the Pn phases by explosion source functions estimated in earlier studies. This leads to azimuthally varying elastic structures that are implemented in the teleseismic intercorrelation procedure. Resulting yield estimates are not significantly changed by the use of the independently estimated layered source structure versus the earlier use of half-space models with delayed pP to account for some overestimate of shallow velocity, largely due to the robustness of the intercorrelation procedure for explosion events in close proximity. The source function deconvolution and layered modeling procedure used here improve on earlier half-space assumptions, but still represent a simplified version of the actual structure and overall source model. Further efforts to include 3D Green's functions for rough near-source topography may increase confidence in the yield estimates but likely will not lead to significantly different results for nearby events.

Data and Resources
Body-wave recordings from global seismic stations were downloaded from the Incorporated Research Institutions for Seismology data management center (IRIS-DMC; http://ds.iris .edu/wilbert3/find_event).

Data Availability Statement and Declaration
All data used in this study are openly available from the sources listed in Data and Resources.

Declaration of Competing Interests
The authors acknowledge that there are no conflicts of interest recorded.
article. Thorne Lay's research on nuclear explosion monitoring is supported by the Air Force Research Laboratory under Contract Number FA9453-18-C-0065.