Characteristics of infrasound signals from North Korean underground nuclear explosions on 2016 January 6 and September 9

Junghyun Park,1 Il-Young Che,2 Brian Stump,1 Chris Hayward,1 Fransiska Dannemann,1,* SeongJu Jeong,1 Kevin Kwong,1 Sarah McComas,1 Harrison R. Oldham,1 Monique M. Scales1,† and Vanshan Wright1 1Roy Huffington Department of Earth Sciences, Southern Methodist University, P.O. Box 750395 Dallas, TX 75275-0395, USA. E-mail: junghyunp@smu.edu 2Korea Institute of Geoscience and Mineral Resources, 124 Gwahak-ro Yuseong-gu, Daejeon 34132, Korea

is supplemented by the addition of regional infrasound arrays such as those in Central Europe (Le Pichon et al. 2008), the Western United States , and the Korean Peninsula (Che et al. 2009).
Direct infrasound from a contained explosion or an earthquake is generated when near-source ground motions at the ground surface are large enough to generate low-frequency, atmospheric acoustic waves. Since these atmospheric waves result from near-source processes, the signals are smaller than those from atmospheric explosions but when observed they can contribute to discrimination and characterization of the source as they are sensitive to source type, depth and size (Arrowsmith et al. 2012;Jones et al. 2015).
North Seismo-acoustic arrays across southern Korea that are jointly operated by the Korea Institute of Geoscience and Mineral Resources (KIGAM) and Southern Methodist University (SMU) recorded seismic and infasound data that can be used to charcterize these events. The explosions were well recorded on seismic elements at the arrays, while the quality and number of infrasound signals are more variable due to atmospheric conditions, source size/depth and noise characteristics at the time of the explosions. The 2006 explosion had no distinct infrasound observations, but the 2009 and 2013 explosions produced clear infrasound signals that were used in studies of infrasound detection, location, propagation modelling as well as yield estimation (Che et al. 2009(Che et al. , 2014. Recent work by Assink et al. (2016) compares infrasound detections and propagation modelling based on atmospheric data for the 2013 and 2016 explosions using two nearby IMS infrasound arrays operated by CTBTO in Japan and Russia.
In this paper, infrasound data recorded at seismo-acoustic arrays in the southern Korean Peninsula are described in Section 2. Signal detections based on automatic detectors are used to identify the infrasound phases and their characteristics as described in Section 3. In Section 4, infrasound propagation models for UNEJ16 and UNES16 are used to assess and interpret the infrasound observations and document the importance of seasonal and diurnal variations in atmospheric conditions. The importance of atmospheric models in improving event location estimates is discussed in Section 5. Finally, in Section 6, the equivalent infrasound energy release for UNEJ16 is estimated using empirical scaling equations for surface explosions and is followed by a summary in Section 7. These analyses illustrate the contributions that infrasound observations can provide in refining the characterization of UNEs.

DATA
KIGAM and SMU cooperatively operate six, seismo-acoustic arrays with varying geometries in South Korea, namely, BRDAR, CHNAR, KSGAR, KMPAR, TJIAR and YPDAR (Fig. 1a). Two additional arrays, ULDAR and YAGAR, are operated by KIGAM. This dense network provides a basis to detect and locate surface explosions at local and regional distances Che et al. 2014). KMPAR has a sampling rate of 100 sample s -1 while the other arrays are sampled at 40 sample s -1 . All acoustic gauges in the arrays are Chaparral Physics Model 2 microphones. These have a flat response from 0.1 to 200 Hz. The signals are recorded using 24-bit digitizers. In order to characterize the wavefield, the acoustic arrays have multiple apertures ranging from 0.2 to 1.0 km with colocated seismometers at some array locations (Fig. 1b). The arrays have co-located weather sensors installed 2 m above the surface that record wind velocity, wind direction and temperature (4 sample s -1 ). Acoustic gauges are each attached to 10, 8-m porous hoses in a star configuration designed to reduce background noise from winds at the boundary layer.

I N F R A S O U N D O B S E RVAT I O N S
Infrasound associated with contained explosions is generated by three mechanisms at different locations: epicentral infrasound generated from surface movement directly above the source; diffracted infrasound generated along the path from the event to the receiver as seismic waves interact with topographic or geological features; and local infrasound generated by seismic ground motion near the receiver . Che et al. (2014) illustrate local, diffracted and epicentral infrasound signals observed at KSGAR from three previous North Korean UNEs in 2006, 2009 and 2013. In the cases of the two UNEs in 2016 discussed in this paper, KS-GAR recorded both local and epicentral infrasound.
Infrasound arrivals associated with UNEJ16 and UNES16 are displayed in Fig. 2 which shows beamformed infrasound waveforms along azimuths consistent with the location of the test site filtered from 1 to 7 Hz at each array. Observations that have identified infrasound signals are included in the figure. Signal detections based on the progressive multichannel correlation (PMCC; Cansi 1995) technique are displayed with 4-element infrasound waveforms (Fig.  3). Automated detection processing parameters using this approach include a 20 s time window, a 50 per cent window overlap, and a 0.1 s consistency threshold. The family settings that group nearestneighbour detections with similar characteristics in the time and frequency domain use group velocities from 0.220 to 0.400 km s -1 . After verifying that the resulting automated infrasound signal detections were consistent with the source direction, they were refined by analyst review using frequency-wavenumber (F-K) analysis in Geotool (Coyne & Henson 1995). Arrival time (first arrival at array elements of detected signal), back-azimuth, correlation, celerity (based on seismic origin time), phase velocity, signal-to-noise ratio (SNR), F-statistic and signal duration for all arrivals identified are summarized in Table 1. Observations at increasing range show decreased signal correlation, SNR, F-statistics and duration, and increased back-azimuth deviations.
Celerity and phase velocity can be used to identify arrivals associated with epicentral infrasound based on previous studies (Che et al. 2009;Negraru et al. 2010;Assink et al. 2012;Che et al. 2014). The range of celerities used to identify tropospheric (Iw), stratospheric (Is) and thermospheric (It) arrivals are >310, 270-310 and <270 m s -1 , respectively, capturing seasonal atmospheric variations. Phase velocity estimates depend on the turning height of the ray (Green et al. 2010) and thus can also help constrain arrival type. The It phases can have a relatively higher phase velocity as a result of the steep incident angle of the wave front at the array. The expected phase velocity for each duct strongly depends on the ducting conditions for the given propagation path. When no The geometry and aperture of each Korean array is plotted. Black circles and grey crosses represent the infrasound sensors and seismometers, respectively. For YAGAR and ULDAR, only the infrasound array information is summarized here with the full-array described by Che et al. (2009).

Figure 2.
The beam-formed infrasound record sections in the direction of the test site (1-7 Hz) from the two North Korean explosions, (a) UNEJ16 and (b) UNES16, where stations are plotted with increasing range. The blue dashed line is a seismic group velocity of 7 km s -1 . The first and second green dashed lines represent the infrasound celerities of 0.4 and 0.2 km s -1 . Signal arrival times estimated by the automatic detection processing (TJIAR result of UNEJ16 from analyst review) are marked as red lines and identified as stratospheric (Is) or thermospheric (It) arrivals. Each waveform in the infrasound celerity range is expanded in Fig. 3. stratospheric duct is present, lower phase velocities can be expected for It phases. Therefore, the It phase may not necessarily have a higher phase velocity than the Is phase, except when considering propagation along the same azimuth as the Is phase.
Based on these criteria, Is and It arrivals from UNEJ16 were identified across the Korean arrays, while only It arrivals were identified from UNES16 (Table 1). The phase velocity estimates for the Is phases span from 323 to 354 m s -1 while those for the It phases range from 329 to 356 m s -1 for both explosions illustrating the phase velocity overlap for these arrivals and the importance of celerity estimates to separate the two. Lower phase velocities can sometimes be expected for the It phase, as documented by Assink et al. (2012).
In the case of UNEJ16, there is a lower phase velocity estimate (323 m s -1 ) for the Is phase compared to the It phase (342 m s -1 ) at KSGAR with the separation of the two phases based on celerity. The closest stations, KSGAR and CHNAR, have similar detection features suggesting common ducting heights. Multiple ducting conditions can produce both stratospheric and thermospheric returns at a particular station (Che et al. 2011;Assink et al. 2012;Fee et al. 2013). The second arrival at CHNAR from UNEJ16 has a reduced celerity (266 m/s) and high amplitude, possibly indicative of a second Is phase. The arrival may be associated with anomalous propagation conditions from a Sudden Stratospheric Warmings (SSW) (Donn & Rind 1971;Assink et al. 2014) that creates a lower mesospheric duct instead of a stratospheric duct and results in an unusually low celerity. Unusual arrivals with abundant high-frequency energy and extremely low celerity (<∼0.19 km s -1 ), have been identified and associated with propagation in the mesosphere (Fee et al. 2013). Based on celerity estimates, KMPAR  Infrasound detections from the underground nuclear explosions (UNEs) in 2016 January 6 (UNEJ16) and 2016 September 9 (UNES16), estimated using the automatic detector PMCC followed by analyst review based on the F-K analysis using Geotool.  Assink et al. (2016). Arrival times and backazimuths used in the location processing are marked by a star next to the Phase ID.
has two identified Is arrivals while the two island arrays (YPDAR and ULDAR) each have a single Is arrival. Arrivals at YAGAR and BRDAR identified as It from UNEJ16 have a lower frequency content (0.1-2.0 Hz), lower SNR and relatively lower phase velocity estimates ( Table 1). The distant station TJIAR (566 km) has an ambiguous automated signal detection that was identified by an analyst as a potential Is phase. Assink et al. (2016) identified one arrival from UNEJ16 at IS45 that can be associated with either a tropospheric or stratospheric return (Table 1), but no observation at IS30. The lack of signal at IS30 might be due to the long distance from the source (1183 km), local noise, and atmospheric conditions at the time of the explosion or a combination of these factors.
In the case of UNES16, the celerity range of infrasound arrivals is relatively low, 259-272 m s -1 (except for the first arrival at YAGAR) with higher phase velocity estimates (347-356 m s -1 ) indicative of It phases (Table 1). The detected signals at KSGAR and YAGAR have a lower frequency content than those from UNEJ16 (Fig. 3). Typically, It phases have lower frequency content than Is phases, due to a mixture of nonlinear propagation as density decreases in the thermosphere and increased attenuation in the thermosphere (Assink et al. 2012;Lonzaga et al. 2015). The It arrivals recorded at these arrays have higher values of signal correlation and SNR and a longer signal duration under lower noise conditions. Some of the phases identified as It from the second event are also observed at higher frequency (1-7 Hz) in Fig. 2(b). The first arrival recorded at YAGAR has a higher celerity (282 m s -1 ) but is also interpreted as an It phase based on its low frequency content of 1-2 Hz. This arrival might be a mesospheric arrival caused by scattering from mesospheric inhomogeneities which are not included in available global atmospheric models (Chunchuzov et al. 2011;Assink et al. 2012).
UNEJ16 epicentral infrasound amplitudes are nearly twice as large as signals from previous explosions that have a similar magnitude and may be related to the atmospheric conditions at the time of the explosion. Seismic arrivals also can couple to acoustic channels (local infrasound) when the local ground motions are large enough that the velocity of microbarometer produces a pressure signal (Kim et al. 2004(Kim et al. , 2010 as observed at KSGAR, YAGAR and CHNAR for P waves with a group velocity close to 7 km s -1 (record sections in Fig. 2). KSGAR, YAGAR and CHNAR, which are closest to the source (307, 350 and 376 km, respectively), have the largest local infrasound signals from both UNEJ16 and UNES16, with relatively high correlation for two distinct arrivals. The peak-to-peak amplitudes of local infrasound (filtered 1-7 Hz) from UNEJ16 are 0.511, 0.423 and 0.285 μbar for KSGAR, YAGAR and CHNAR, respectively, while those from UNES16 are 0.643, 0.211 and 0.350 μbar.

I N F R A S O U N D P RO PA G AT I O N M O D E L L I N G
Infrasound propagation is dependent on temperature, wind speed and direction. The combined effects of wind and temperature can be represented by the effective sound speed approximation (Godin Downloaded from https://academic.oup.com/gji/article-abstract/214/3/1865/5042944 by Southern Methodist University Libraries user on 12 July 2020 2002). Under typical conditions the effective sound speed approximation is appropriate for assessing signals that return to the Earth from the stratosphere and below (Assink et al. 2017). Realistic atmospheric data and models provide input parameters necessary to calculate the effective sound speed for a stratified atmosphere, where γ is the ratio of specific heats, R is the gas constant for air, T is the absolute temperature andn · u projects the wind ( u) in the direction from source to observern (Negraru et al. 2010).
Ground-to-Space (G2S) atmospheric specifications (Drob et al. 2003) at or near the time of the explosion are used to compute the effective sound speed. These models are based on the well resolved and constrained operational meteorological analysis fields from the National Oceanic and Atmospheric Administration (NOAA) Global Forecast System (GFS) analysis fields below 35 km (Kalnay et al. 1990), the National Aeronautics and Space Administration (NASA) Modern-Era Retrospective analysis for Research and Applications analysis fields (Rienecker et al. 2011) between 25 and 75 km, as well as the less well constrained National Research Laboratory (NRL) Mass Spectrometer Incoherent Scatter (MSIS R )/Horizontal Wind Model (HWM14) empirical upper atmospheric model at higher altitudes Drob et al. 2015). The hourly zonal and meridional winds and associated effective sound speeds for 00-23 hr local time on the day and location where both explosions were conducted (OT: 01:30:01 UTC 2016 January 6, UNEJ16 and 00:30:01 UTC 2016 September 9, UNES16) are plotted in Fig. 4. The 24, hourly profiles at the same location (grey lines) for the days of UNEJ16 and UNES16 are compared with the averaged profile (black lines) and the profile for the time of two explosions (red lines) in Figs 4(a), (b), (e) and (f). The hourly atmospheric profiles below an altitude of 50 km are similar, while wind speeds above this altitude document hourly variations as a result of solar heating. Note that the upper part of G2S specification has an uncertainty associated with the semi-empirical MSIS and HWM models which are based on a long dataset of many years of observations but does not rely on actual assimilated observations. The effective sound speeds for all arrays were calculated using the atmospheric profiles at 01 and 00 UTC, close to the respective times when UNEJ16 and UNES16 were detonated (Figs 4c and g).
To explore infrasound detectability as a function of azimuth and distance under the assumption that the atmosphere is horizontally stratified and independent of range, the effective sound speed ratios (c eff −ratio ; ratio of effective sound speed to sound speed at the ground) for both explosions are plotted as a function of altitude (Figs 4d and h). Sound ducting is predicted (favourable downwind propagation), where c eff −ratio is greater than one. In the case of UNEJ16, the paths to the two IMS infrasound arrays (IS30 and IS45) show favourable conditions for eastward tropospheric arrivals, while all Korean arrays are in a region that is slightly favourable for stratospheric propagation turning at approximately 50 km altitude. For the UNES16, both tropospheric and stratospheric propagation is unfavourable to all the Korean arrays and IMS stations but thermospheric ducting above 100 km altitude may explain the detections.
Following the method of Blom & Waxler (2017), 3-D ray tracing was conducted for the time nearest the detonation time for UNEJ16 and UNES16 (Fig. 5). The ray tracing calculation goes beyond the effective sound speed approximation, and rigorously treats the effects of temperature and wind within the limitations of the ray approximation. We use a single G2S profile at the source, with a source elevation of 2.0 km based on the average source and array elevations. This method approximates geometric ray paths and estimates relative amplitudes resulting from an impulsive signal in an inhomogeneous moving medium (Blom & Waxler 2017). Ray tracing using spherical coordinates was completed for inclination angles from 0.5 • to 60.5 • relative to the horizontal at a step of 1 • , azimuths from 0 • to 360 • with a step of 2 • and a maximum of 5 surface bounces. A frequency band of 1 Hz and an absorption coefficient of 0.3 were initially used to focus on Is arrivals. The details of the ray tracing parameters are described in the GeoAc manual (Blom 2014).
The predicted ray turning heights are plotted in Fig. 5 for tropospheric (altitude < 15 km), stratospheric (15 km < altitude < 80 km), and thermospheric (80 km < altitude) arrivals. The majority of the predicted arrivals from both explosions are thermospheric across all azimuths, although the amplitudes are expected to be low as a result of attenuation. In the case of UNEJ16, tropospheric arrivals are predicted NE to S of the source, with Is arrival propagation to the S, N and W. The Korean infrasound arrays that are to the SW and S of the explosion are in a region with predicted It arrivals as well as Is arrivals but no predicted Iw arrivals ( Fig. 5a). Iw, Is and It phases are predicted at IS45, while no arrivals except for possible low amplitude multiple It arrivals at IS30 were predicted. These results are consistent with the study of Assink et al. (2016) who found that the propagation paths for the Iw and Is phases have similar travel times, which is unusual, as these are typically separated in time. There are bi-directional stratospheric ducting conditions, likely due to the occurrence of a minor SSW (Assink et al. 2014). It was also found that the number of detections from the E at CHNAR were reduced during the winter due to a mixture of eastward and westward propagation from a SSW (Park et al. 2016). In the case of UNES16 (Fig. 5b), Is arrivals are predicted to the W and NW, directions that are unfavourable to both the Korean and IMS arrays. Even though thermospheric predictions occur at all azimuths, only the closest three stations, KSGAR, YAGAR and CHNAR, have possible thermospheric returns. Celerities estimated from the data at these stations are relatively low (Table 1), consistent with It arrivals. The similarities and differences between the observations and the predictions will be discussed in more detail along with eigenray results from the ray tracing calculations in Section 5.
The relative amplitude, defined as the ratio of the amplitude in the far field relative to that at 1 km from the source assuming frequency-independent geometric spreading, is calculated based on the attenuation model of Sutherland & Bass (2004). The relative amplitudes along the ray paths decrease with distance from the source (Fig. 6). Relative amplitudes from -200 to 0 dB are plotted in the figure and the average relative amplitudes for Iw and Is phases are summarized. Note that the Eikonal equation from linear theory used in this study provides estimates of arrival time even in the nonlinear regime but the amplitudes must be estimated separately as outlined in Lonzaga et al. (2015). As these authors note, infrasound attenuation in the thermosphere requires a consideration of nonlinear propagation conditions, which influence both the observed frequency and acoustic pressure of the arrival (e.g. Lonzaga et al. 2015). The linear calculations reported here do not include these effects and thus may dramatically overestimate the wavefield attenuation, thus the It arrival predictions are only used as a comparative guide to the observations.
In the case of UNEJ16, strong relative amplitudes for Is arrivals are predicted to the NW and N at regional distances, where there are no stations. Tropospheric energy is dominant to the NE (up to 1200 km), E (within 1000 km), SE and S (within 400 km). Relatively high average relative amplitudes for It arrivals are predicted to the SW and S, where the Korean infrasound arrays are located (i.e. KSGAR). IS30 has a predicted It arrival with a low relative amplitude, consistent with the lack of observations at the array. IS45 has relatively high energy for both tropospheric and stratospheric returns but lower than for the Korean arrays. In the case of UNES16, no tropospheric energy was found and stratospheric returns are dominant to the NW (-71.53 dB) and W (-73.21 dB), with a slightly reduced (∼2 dB) average relative amplitude prediction compared to UNEJ16. For this second explosion, the energy level for It arrivals are predicted to be a relatively high to the N, NW, W and SW. The Korean infrasound arrays are located in a region of stronger predicted amplitudes for the thermospheric returns in the case of UNEJ16, while these are on the edge of the ensonified region for UNES16. However, only the three closest stations from the source (KSGAR, YAGAR and CHNAR) have observed arrivals which may reflect the unaccounted nonlinear effects in the thermosphere previously discussed. The two IMS stations have predicted It arrivals that are lower energy for this event.
In order to explore the expected year-to-year variability of atmospheric models for the time periods of UNEJ16 and UNES16, Downloaded from https://academic.oup.com/gji/article-abstract/214/3/1865/5042944 by Southern Methodist University Libraries user on 12 July 2020  Downloaded from https://academic.oup.com/gji/article-abstract/214/3/1865/5042944 by Southern Methodist University Libraries user on 12 July 2020 ray tracing was conducted for all possible G2S atmospheric models from 2006 through 2016 for January and September. Atmospheric specifications at 00:00:00 UTC were used because older G2S datasets are only available at 6-hr intervals. These simulations used the same ray tracing parameters as before, with a maximum number of surface bounces of 5 and a dominant frequency of 1 Hz. The highest relative amplitudes were extracted with a cut-off value of -25 dB. Fig. 7 displays the density distributions for the number of ray tracing hits for each month [(a) January and (b) September] based on ray tracing through all the models covering the 11 yrs. The resulting distributions are compared to the day specific model predictions for (a) UNEJ16 and (b) UNES16 (colour plots) from Fig. 6. In January, there are dominant tropospheric returns to the NE, E, SE and S and stratospheric returns to the S-W-N at the time of UNEJ16 (Fig.  7a). Due to the early stage of a SSW in January 2016, stratospheric returns are unusually dominant to the west, which typically occurs during the summer. Based on the historical ray tracing results, eastward stratospheric predictions are significant during this month. These models indicate that abnormal atmospheric conditions during the January 2016 provides favourable conditions for detection at the infrasound stations on the Korean peninsula. Thermospheric returns are predicted to the SW, where KMPAR, YPDAR and BRDAR are located. The observations from the explosion are consistent with the ray tracing density distribution for January 2016 and illustrate the predictive power of the historic atmospheric models. Note that the thermospheric distributions have relatively minor variations, as these ray paths are sensitive to the semi-empirical HWM/MSIS specifications that are more or less the same each year except for variations due to space weather (Garcés et al. 2002). In addition, stratospheric energy for this month is enriched to the E and SE, where there are no observations at the time of the explosion. Comparisons suggest that strong tropospheric returns are expected in the region covering from the NE to S (clockwise), while stratospheric returns are unusually favourable in all directions during January. This observation is consistent with the study by Park et al. (2016) that documented the temporal detections from both E and W due to a SSW during the winter in this area. Stratospheric returns which are more dominant than thermospheric returns were unfavourable to the E and SE at the specific time of UNEJ16. Predicted thermospheric returns for UNES16 have higher predicted energy and are consistent with the high signal correlation values, F-statistics and SNR, and long signal duration (Table 1) from detections at the closest Korean infrasound arrays (Fig. 7b). Propagation characteristics in September are more unidirectional than in January as this month is not characterized by the occurrence of SSWs. Single predictions as well as the density estimates indicate that neither tropospheric nor stratospheric paths to the Korean infrasound arrays or the IMS stations are expected, as a result of the equinox period. Tropospheric predictions towards the Korean peninsula are unlikely for either explosion, as these are mainly controlled by the direction of the jet stream which is primarily eastward, especially during January (Park et al. 2016).
In order to further investigate the temporal variations in the atmosphere, ray paths including turning heights and relative amplitudes were computed using G2S profiles for every hour (0-23 hr, local time) on 2016 January 6 and September 9. Fig. 8 shows the hourly ray tracing that include the Korean arrays (SW and S), IS45 array (NE) and IS30 array (SE) for both UNEJ16 and UNES16. Ray tracing predictions vary in both space (station location) and time (hourly). Both explosions in 2016 were conducted in the morning (10:30:00 am for UNEJ16 and 09:30:00 am for UNES16).
In the case of UNEJ16, KMPAR and YPDAR have predicted It arrivals during the entire day, while this arrival is only predicted during the morning at CHNAR and BRDAR (Fig. 8a). Is arrivals are predicted only to the S, including KSGAR, YAGAR and TJIAR, during the morning. ULDAR has a lower expectation for this arrival (Fig. 8b), but could have a tropospheric arrival if the explosion was detonated in the afternoon. IS45 has both tropospheric and Is arrivals during the entire day (Fig. 8c) while IS30 has a predicted Is arrival during the night (Fig. 8d).
For all hours of the day, UNES16 has only thermospheric predictions (bottom figures of Fig. 8). CHNAR, KSGAR and YAGAR have thermospheric predictions only during the morning (Figs 8e and f). KMPAR, YPDAR and BRDAR to the SW have predicted It arrivals at the explosion time, but there were no observations (Fig.  8e). These arrays have predicted It arrivals during the late night. It arrivals are predicted at ULDAR and IS30 during the early morning and night, and at TJIAR and IS45 during the early morning (Figs  8f-h). These hourly variations in the ray tracing predictions illustrate the importance of time varying atmospheric conditions over the course of a day. Che et al. (2018) demonstrate that detection capabilities in this region are relatively stable during summer due to steady stratospheric wind conditions, with large variations from September to May when stratospheric winds are more variable. Both explosions were conducted during the time period (January and September) when the atmospheric structure may be more variable, suggesting that the detectability of Is arrivals at the stations might also be expected to be variable in time. The hourly variation shown in Fig. 8 is mostly due to the atmospheric dynamics in the mesosphere and lower thermosphere (e.g. atmospheric tides) (Le Pichon et al. 2005;Assink et al. 2012).
These hourly variations in the ray paths are provided as supplementary movie files. The models document significant changes in azimuthal patterns of It arrivals over the 24-hr period for the both the January and September explosions. In the case of UNEJ16, the predicted pattern of tropospheric arrivals rotates counter-clockwise by about ∼10 • over the 24-hr period, while those for stratospheric rays show a larger clockwise rotation of ∼30 • . The model predictions for UNES16 support stratospheric returns over a small area to the west, from 6 am to 12 pm. The detectability of infrasound signals along the edges of the illuminated zones may be sensitive to the precise time of the source as illustrated by these time dependent rotations.

I N F R A S O U N D E V E N T L O C AT I O N
The modelling results illustrate the strong degree to which infrasound propagation is dependent on temporal variations in the atmosphere. Accounting for the associated variations in travel times and back-azimuths offers the opportunity to improve infrasound event location estimates (Blom et al. 2015). Eigenrays were identified from the ray tracing at the time nearest the explosion for each station in order to estimate azimuthal corrections for application to the location process. Estimated eigenrays for UNEJ16 and UNES16 are plotted in Fig. 9 with details summarized in Table 2.
Comparing the observations (Table 1) with the model predictions (Table 2), there are some similarities and differences as summarized below.
(1) There are more predicted It eigenrays for the thermospheric returns than observations for both explosions. For the UNEJ16 case, there are predicted It arrivals at all stations except for BR-DAR and YPDAR. BRDAR has no eigenray, but a possible It  arrival was identified. In the case of UNES16, all stations have eigenrays (Table 2), while the It arrivals were only recorded at KS-GAR, CHNAR and YAGAR (Table 1). These differences might reflect the linear nature of the amplitude calculations as noted earlier.
(2) For both explosions, the three closest arrays, KSGAR, YA-GAR and CHNAR, have It phase detections with celerities higher that the predictions. This comparison suggests that It arrivals at these stations may have lower ducting altitudes than the predictions. These arrivals might be a mesospheric arrival generated by sound scattering from mesospheric inhomogeneities which is not captured by the current global atmospheric model (Chunchuzov et al. 2011;Assink et al. 2012).
(3) In the case of UNEJ16, Is phases were observed at six infrasound arrays. There were no Is arrivals from UNES16. The observations at KSGAR, YPDAR and TJIAR, are well matched with Is predictions, while no predictions were found for CHNAR, KMPAR and ULDAR. Lonzaga et al. (2015) demonstrate that modelling using currently available G2S atmospheric specifications can result in no predicted arrivals due to uncertainties in the profile as a function of increasing altitude and argue that these models may need small-scale updates to the atmospheric profile.
(4) One observation associated with UNEJ16 with a celerity of 322 m s -1 was found at IS45 (Assink et al. 2016), that matches four predicted eigenrays (one for the stratospheric, and three for the tropospheric returns with a maximum of 5 bounces). Ray tracing produces one additional thermospheric prediction, but this phase was not observed. Multiple eigenrays for the It phase were predicted at IS30, while no observations were found, suggesting that the energy was sufficiently attenuated to obscure the arrival.
Infrasound observations (Table 1) consistent with the predictions were selected for location processing with observed azimuths corrected by the maximum back-azimuth deviation from the eigenray calculations. Arrival times and back-azimuths used in this processing are marked by stars in Table 1. The associated eigenray information used for azimuth correction are also marked by stars in Table  2. Detections at IS45 for UNEJ16 from Assink et al. (2016) were included in the location process. Event locations are estimated using the Bayesian Infrasonic Source Location (BISL; Modrak et al. 2010) methodology, accounting for unknown source-to-array path effects with a uniform probability for infrasound group velocity. In this study, group velocities from 0.24 to 0.35 km s -1 were used. The assumed standard deviations for azimuth and arrival time were bounded between 1 • and 3 • , and 10 and 30 s in order to bound the error ellipses.
Event locations without wind corrections to the measured azimuths are compared to locations with wind corrected azimuths based on the ray tracing simulation in Fig. 10, and include the 95 per cent credibility contours for each location. The seismic origin times and locations provided by the USGS are compared to the infrasound estimates from BISL for both explosions (Table 3). The UNEJ16 location without wind correction is biased to the east of the seismic location (up to 28 km), indicating that initial observations are affected by tropospheric and stratospheric winds. The initial UNES16 location is biased to the north of the seismic location (up to 165 km), because of the limited number of detections consistent with atmospheric winds at the time of the explosions and the poor station coverage to the north and east of the explosion.
In the case of UNEJ16, the distance between the seismic and infrasound location estimates is reduced by 50 per cent with the wind corrected azimuths (from 26 to 11km) using the larger standard deviations for azimuth and arrival time (3 • and 30 s). Similar locations for UNES16 produce distances between the seismic location and the BISL estimates with much larger (up to 5 times) uncertainty due to the fewer number of arrivals. For the UNEJ16, the location estimate using the larger standard deviations of azimuth and arrival time (3 • and 30 s) provides a location that is closer to the seismic origin time and surface location than the infrasound location calculated with smaller standard deviations. This result suggests that under certain conditions the resulting location estimate can be biased in order to reduce the uncertainty contour. After wind correction, estimates of origin time and location for UNES16 are closer to the seismic estimates with the smaller standard deviations (1 • and 10 s).
These results illustrate the utility of including azimuth corrections based upon the atmospheric wind profile at the time of the explosion as well as the need to improve estimates of standard deviation estimates for back-azimuth and arrival time. It may be possible to modify travel time priors as well based on the time dependent atmosphere in order to further decrease the estimated error ellipses (Blom et al. 2015). Mutschlecner et al. (1999) developed an empirical relation between wind conditions and infrasound amplitudes based on observations from Nevada Test Site nuclear tests (1951)(1952)(1953)(1954)(1955)(1956)(1957)(1958) that provides a methodology to estimate the yield of the equivalent atmospheric explosion. They argue that accounting for stratospheric winds is crucial to relating the observed infrasonic amplitude to this effective explosion yield. Whitaker et al. (2003) applied the empirical scaling to yield estimates of equivalent high explosive (HE) nearsurface explosions. The relation has been applied to a number of explosions including previous North Korean UNEs (Che et al. 2009(Che et al. , 2014, the Buncefield Oil Depot explosion (Ceranna et al. 2009), and a gas-pipeline explosion (Evers et al. 2007). Even though the infrasound energy in the case of a contained explosion is generated by near-field strong motion at the free surface above the underground explosion, this assessment provides a comparison of the size of the equivalent explosive yield from the associated infrasound arrivals. This equivalent yield estimate is related to the absolute size of the contained explosion, its depth of burial, local material properties, and topography. As a result, the yield scaling for contained and surface explosions might be expected to be different and a subject for future research. Table 2. Predicted infrasound arrival characteristics at each array for underground nuclear explosions (UNEs) at 01:00:00 UTC 2016 January 6 (UNEJ16) and at 00:00:00 UTC 2016 September 9 (UNES16) based on ray tracing. Note: The associated eigenray information used for azimuth correction are marked by a star next to the Phase ID.

I N F R A S O U N D S O U RC E E N E RG Y
In order to estimate the equivalent atmospheric explosion yield associated with the UNE, the peak-to-peak pressures of Is arrivals identified at CHNAR, KSGAR, KMPAR, ULDAR, YPDAR and TJIAR for UNEJ16 were measured. The two separate Is arrivals at both CHNAR and KMPAR were included in the calculation. Data from UNES16 could not be used to estimate infrasonic source energy using this methodology, since only It arrivals were observed. The stratospheric wind corrected pressure amplitude (P wca , μbar) is defined as, where P raw is the raw peak-to-peak pressure amplitude of the stratospheric signal, k is an empirical constant (k = 0.0018, s m -1 ), and V d is the wind component directed from the source to the array based on the G2S specification at altitude of 50 km. The wind speeds were taken at the mid-point between source and receiver (Table 4). This empirical equation was developed using a least squares regression of wind corrected pressure amplitudes and scaled range (Whitaker et al. 2003): P wca = 59457 * (S R) −1.4072 , where S R is the scaled range (S R = R/(2 × chgwt) 0.5 ), R is the source to receiver range in kilometres and chgwt is the charge weight in metric kt of HE. The wind corrected amplitude from the stratospheric observations at the six arrays range from 0.1605 to 1.0431 μbar. Individual measurements were used to estimate source strengths from 0.9 to 16.1 tons of equivalent TNT with an average value of 6.4 tons and standard deviation of 4.6 tons (Table 4). ULDAR has a significantly compensated value after the wind correction because the stratospheric wind is opposite to the site (-4.0 m s -1 ) at altitude of 50 km, and produces the highest yield, 16.1 ton. On the other hand, KMPAR, at a similar range to ULDAR, has reduced amplitude due to a strong stratospheric wind speed (9.2 m s -1 ). These results illustrate that this yield estimation is sensitive to the stratospheric wind speed and direction, although the eigenrays of  (UNEJ16 and UNES16), is marked as a red star. The geometry and aperture of each Korean array are described in Fig. 1(b). Each figure shows the average infrasound event location for UNEJ16 with standard deviations of (a) 3 • and 30 s and (b) 1 • and 10 s and for UNES16 with standard deviations of (c) 3 • and 30 s and (d) 1 • and 10 s using the Bayesian Infrasonic Source Location method. Locations with uncorrected (grey circle) and corrected (blue circle) azimuths based on a ray tracing are plotted with associated 95 per cent confidence contours. An expanded source location area is shown in the left corner of each figure. The uncorrected and corrected azimuths are plotted from arrays as grey and blue lines, respectively. The detection back-azimuth at IS45 for UNEJ16 is from Assink et al. (2016). Table 3. The seismic origin times and locations from the USGS for underground nuclear explosions (UNEs) on 2016 January 6 (UNEJ16) and 2016 September 9 (UNES16) are compared to estimates of event locations and origin times from Bayesian Infrasonic Source Location (BISL; Modrak et al. 2010)  Notes: The BISL calculations with and without wind corrections were made with two sets of standard deviations (SD) for azimuth and arrival time. Distances between the USGS seismic locations and the infrasound locations are summarized in the 'DIS' column. Table 4. Yield estimates for underground nuclear explosion (UNE) on 2016 January 6 (UNEJ16), based on the peak-to-peak pressures (P raw ; raw pressure amplitude, and P wca ; wind corrected pressure amplitude) of stratospheric arrivals corrected for stratospheric wind speed (V d ) at 50 km altitude using the Ground-to-Space (G2S) specification (Drob et al. 2003 KSGAR, YPDAR and TJIAR have turning heights of ∼50 km (Fig. 9). Taking the stratospheric wind at an altitude 50 km may not be appropriate in this case as the stratospheric duct is likely atypical from the normal summer situation in which 50 km is a reasonable estimate for the top of the stratospheric waveguide. Moreover, the stratospheric winds at this altitude during the winter are far more uncertain than during the summer, especially since it is during the onset of a SSW (Assink et al. 2016). Le Pichon et al. (2015 demonstrated that the uncertainty of wind profiles increases with altitude above 40 km by comparing the observations with several model profiles. Fig. 11 shows the sensitivity of the yield estimation to stratospheric winds in a given range of refraction altitudes (40-60 km) with their associated uncertainties. The wind profiles are taken from the mid-point between source and receivers. We also find that the profiles tend to be highly variable in the altitude range between 48 and 52 km (Fig. 11a). Unfavorable wind condition to all stations (around -10 m s -1 ) below 48 km greatly impacts the wind corrected pressure amplitude (average of ∼1 μbar) and the resulting yield estimates (average of 15 ton) (Figs 11b and c). UL-DAR has the highest yield estimate (up to 32 ton) using the lowest wind speed at 47 km. On the other hand, above 52 km, consistently strong wind to all stations provides reduced pressure estimates (average of ∼0.4 μbar) with an average yield of 3 ton. These results illustrate that yield estimation from infrasound data is strongly dependent on the wind conditions including refraction altitude and wind speed and these parameters may vary between different atmospheric specifications, that is G2S or ECMWF (European Centre for Middle-range Weather Forecast, www.ecmwf.int). The seismic yield of the UNEJ16 is <10 kiloton equivalent yield of TNT (NOR-SAR article, titled in "New nuclear test by North Korea" posteon 2016 May 05, https://www.norsar.no/in-focus/new-nuclear-test-bynorth-korea-article186-863.html) for comparison. Further studies are warranted to not only investigate seismic-to-acoustic coupling as well as the importance of the atmospheric model on the estimation of equivalent yields from infrasound associated with UNEs.

C O N C L U S I O N S
North Korea conducted two underground nuclear tests at the Punggye-ri Nuclear Test Site on January 6 at 01:30:01 UTC and September 9 at 00:30:01 UTC 2016. Infrasound signals from the underground nuclear explosion (UNE) in January 2016 (UNEJ16) are well recorded at eight seismo-acoustic arrays, operated by SMU and KIGAM, across southern Korea, while signals from the second UNE in September 2016 (UNES16) are only recorded at three arrays. Signals from both explosions are used to identify infrasound phases, quantify propagation path effects and estimate the equivalent source strength (only for UNEJ16).
Infrasound signals were identified using PMCC followed by analyst review. Detection results provide estimates of celerity, phase Figure 11. (a) Wind profiles taken at the mid-point between source and receivers; (b) wind corrected peak pressure amplitude, P wca ; and (c) yield estimates as a function of refraction altitude (40-60 km). Two stratospheric arrivals at CHNAR and KMPAR were included in the calculation. The dashed black line is at an altitude of 50 km represents the average profile for estimating P wca and yield. velocity, correlation, back-azimuth, SNR, F-statistic and signal duration at all arrays. Infrasound phase identification was based on celerity since the origin time of the source is well constrained by the seismic data. Phase velocity estimates were also used to verify the phase identification. The Korean infrasound arrays recorded both Is and It arrivals from UNEJ16. Infrasound detections at IS45 are possibly Is or Iw arrivals with no arrivals observed at IS30 (Assink et al. 2016). There were only It arrivals at the closest three arrays from UNES16.
3-D ray tracing was completed using the method of Blom & Waxler (2017) based upon 1-hr atmospheric updates from G2S atmospheric specifications (Drob et al. 2003). Model predictions in terms of range and celerity are compared with observations for both UNEJ16 and UNES16. In the case of UNEJ16, stratospheric winds are marginally adequate for arrivals at the Korean arrays, especially KSGAR and CHNAR. Due to the early SSW in January 2016, stratospheric returns are strong to the west, a condition more typical of summer. The single model prediction nearest the time of the explosion are consistent with average ray tracing predictions from 2006 to 2016. The long-term results also predict strong observations to the east during this month. Note that the upper part of G2S atmospheric specification will be self-similar through the years, because of the use of semi-empirical models in this part of the atmosphere. These broader results indicate that on average in January, there are stratospheric conditions that can contribute to infrasound detection at the Korean arrays. Model predictions for UNES16 are dominated by thermospheric returns at all directions from the source with relatively weak stratospheric returns to the NW. This single G2S model is also consistent with the long-term ray tracing density distributions for September when thermospheric returns are favourable to the Korean arrays, reflecting self-similarity for the thermospheric predictions during September. Che et al. (2018) demonstrate that steady stratospheric winds in this region from June to August provide improved event detectability at the Korean arrays compared to decreased detectability from September to May as a result of highly variable stratospheric winds. Both explosions were conducted during the time periods (January and September) when the atmospheric structure is more variable, suggesting that the detectability for Is arrivals at the stations might be expected to be variable in time. Hourly variations in ray tracing predictions illustrates how atmospheric propagation directions can change over short time periods across the network. For both explosions, the models suggest that along directions where there were no instruments (to the E or NE) that infrasound signals might have been expected. Limited observations at two IMS infrasound arrays (Assink et al. 2016) are consistent with these predictions.
There are similarities and differences between the observations and the eigenray results from the ray tracing. For both explosions, the three close arrays, KSGAR, YAGAR and CHNAR, have It observations with higher celerities than predictions, suggesting that It arrivals at these stations might have lower ducting altitudes than the predictions. These arrivals might be a mesospheric arrivals caused by sound scattering from mesospheric inhomogeneities which are not captured by currently available global atmospheric models (Chunchuzov et al. 2011;Assink et al. 2012). Is phases were observed at six infrasound arrays for UNEJ16. The observations at KSGAR, YPDAR, and TJIAR, are well matched with Is predictions, while no predictions were found at CHNAR, KMPAR and ULDAR. Lonzaga et al. (2015) demonstrated that modelling using currently available G2S atmospheric specification can predict no arrivals as a result of increasing uncertainties in the profile as a function of altitude and suggest there is a need to update small-scale structures in the atmospheric profile in order to capture the missing infrasound arrivals.
Downloaded from https://academic.oup.com/gji/article-abstract/214/3/1865/5042944 by Southern Methodist University Libraries user on 12 July 2020 Estimated event locations using the Bayesian Infrasonic Source Location (BISL; Modrak et al. 2010) methodology are improved by correcting back azimuth measurements based on the ray tracing. Infrasound observations consistent with the predictions were selected for location processing. For UNEJ16 there is a 11-km difference from the seismic epicentre, while the seismic location is included within the formal 95 per cent credibility bounds of the infrasound location. Due to limited detections as a result of atmospheric conditions at the time of UNES16, the initial location is biased by up to 165 km north of the seismic location. After wind corrections to the azimuth estimates, the distance from the seismic epicentre is reduced to 67 km. It may be possible to better constrain the travel time priors using a time dependent atmosphere, and therefore improve the location estimates and further decrease the estimated error ellipses (Blom et al. 2015).
Using an empirical methodology for estimating the equivalent HE for a surface explosion from Is arrivals (Whitaker et al. 2003), infrasound energy release was estimated using all the Is arrivals from UNEJ16. The wind corrected amplitudes for the stratospheric observations at the six arrays ranged from 0.1605 to 1.0431 μbar. The equivalent source energy estimates ranged from 0.9 to 16.1 tons of TNT with an average of 6.4 tons and standard deviation of 4.6 tons. The yield estimates from the infrasound data are sensitive to stratospheric wind speed and refraction altitude.
Whitaker's empirical relation has been applied to previous NK nuclear explosions. For the 2009 explosion, Che et al. (2009) estimated 3.0 (0.7-8.6) tons from observations at KSGAR, YAGAR, ULDAR and BRDAR, using HWM93 (Horizontal Wind Model) (Drob et al. 2008). Based on the amplitude at KSGAR, Che et al. (2014) estimated the equivalent yields of 0.66 tons for the 2009 explosion and 3.39 tons for the 2013 explosion, using the ECMWF atmospheric model. The associated wind-corrected amplitudes for each explosion were 0.18 and 0.57 μbar, respectively. Che et al. (2014) corrected the raw amplitude using stratospheric winds from 25 to 35 km altitude (empirical relation was developed using stratospheric wind at 50 km), based on the predicted refraction height for the Is phase from ray tracing. It should be noted that the atmospheric model has increased uncertainty above 40 km altitude (Le Pichon et al. 2015) and stratospheric winds at this altitude during the winter are far more uncertain than during the summer, especially with the onset of a SSW (Assink et al. 2016). Yield estimates are sensitive to stratospheric winds and refraction altitudes, with profiles that tend to be highly variable. In the case of UNEJ16 below 48 km, there are unfavourable wind conditions to all stations (around -10 m s -1 ) that results in a large wind correction (average of ∼1 μbar) and an average yield estimate of 15 tons. The highest yield estimate is up to 32 tons as a result of the use of the lowest wind speed at 47 km (at ULDAR). Above 52 km, consistently strong wind to all stations leads to a reduced pressure amplitude (average of ∼0.4 μbar) and an average yield of 3 tons.
Direct comparisons of the estimated infrasonic yields for the four North Korean nuclear explosions is hampered by the use of different atmospheric wind models (HWM93/ECMWF/G2S), arrays (or elements), wind-correction calculations and environmental factors. To compare the equivalent surface yield among the four explosions, more work is needed to correct amplitudes at each array using the same atmospheric model(s) and a consistent empirical relation. For example, the raw infrasound amplitudes at KSGAR for the 2009, 2013 and January 2016 explosions are 0.19, 0.68 and 1.27 μbar, respectively (with no observation from the 2006 (Che et al. 2014)). Even though the magnitude of the UNEJ16 (M5.1) is the same as that for the 2013 explosion (M5.1) (USGS 2016 significant earthquakes catalog), the infrasound amplitude for 2016 is almost twice as large as that from 2013. During February 2013, there was a more typical winter situation with predominant eastward stratospheric ducting. There were unfavourable wind conditions to all Korean infrasound arrays, except for KSGAR, at the edge of the stratospheric waveguide, while the more efficient propagation conditions are towards to IS30 and IS45 (Che et al. 2014). In contrast, stratospheric propagation to the SW became more efficient at KSGAR due to the onset of a SSW during January 2016. Full waveform modelling techniques are needed to further explore these effects. In addition, UNES16 cannot be compared with the other explosions since there were no stratospheric observations, even though this explosion had a large magnitude (M5.3). Based on the G2S atmospheric specifications for 2009 and 2013 (4.3 and -28.1 m s -1 at 50 km altitude), the wind corrected amplitudes are 0.15 and 2.18 μbar for the 2009 and 2013 explosions, respectively with infrasonic energies estimated to be 0.6 and 23.5 tons of TNT. This result illustrates the strong stratospheric wind effect in 2013 was accompanied by a significant amplitude correction, and motivates the need to quantify these effects for any relative comparisons of the explosions. Infrasound explosion energy estimates are not only related to the yield of the nuclear test, but also source depth, secondary source effects such as spallation, topography in addition to the atmospheric conditions at the time of the explosion. Finally, the results of this comparative study suggest that detection capabilities for the Korean arrays and nearby IMS stations are seasonally dependent on stratospheric and possibly tropospheric wind dynamics (Che et al. 2012).