A multi-technology analysis of the 2017 North Korean nuclear test

On 3 September 2017 official channels of the Democratic People's Republic of Korea announced the successful test of a thermonuclear device. Only seconds to minutes after the alleged nuclear explosion at the Punggye-ri nuclear test site in the mountainous region in the country's northeast at 03:30:02 (UTC), hundreds of seismic stations distributed all around the globe picked up strong and distinct signals associated with an explosion. Different seismological agencies reported body wave magnitudes of well above 6.0, consequently estimating the explosive yield of the device on the order of hundreds of kT TNT equivalent. The 2017 event can therefore be assessed as being multiple times larger in energy than the two preceding North Korean events in January and September 2016.

This study provides a multi-technology analysis of the 2017 North Korean event and its aftermath using a wide array of geophysical methods. Seismological investigations locate the event within the test site at a depth of approximately 0.6 km below the surface. The radiation and generation of P- and S-wave energy in the source region are significantly influenced by the topography of the Mt. Mantap massif. Inversions for the full moment tensor of the main event reveal a dominant isotropic component accompanied by significant amounts of double couple and compensated linear vector dipole terms, confirming the explosive character of the event. The analysis of the source mechanism of an aftershock that occurred around 8 min after the test in the direct vicinity suggest a cavity collapse. Measurements at seismic stations of the International Monitoring System result in a body wave magnitude of 6.2, which translates to an yield estimate of around 400 kT TNT equivalent. The explosive yield is possibly overestimated, since topography and depth phases both tend to enhance the peak amplitudes of teleseismic P waves. Interferometric synthetic aperture radar analysis using data from the ALOS-2 satellite reveal strong surface deformations in the epicenter region. Additional multispectral optical data from the Pleiades satellite show clear landslide activity at the test site. The strong surface deformations generated large acoustic pressure peaks, which were observed as infrasound signals with distinctive waveforms even at distances of 401 km. In the aftermath of the 2017 event, atmospheric traces of the fission product ¹³³Xe were detected at various locations in the wider region. While for ¹³³Xe measurements in September 2017, the Punggye-ri test site is disfavored as a source by means of atmospheric transport modeling, detections in October 2017 at the International Monitoring System station RN58 in Russia indicate a potential delayed leakage of ¹³³Xe at the test site from the 2017 North Korean nuclear test.

1 Introduction

The Comprehensive Nuclear-Test-Ban Treaty (CTBT) and its associated entity, the Preparatory Commission for the CTBT organization (CTBTO), are dedicated to monitoring and banning nuclear explosions worldwide – underground, in water, or in the atmosphere. The CTBT was opened for signature in 1996 but will only enter into force after the 44 nuclear technology holders (states listed in Annex 2 of the CTBT) will have signed and ratified the treaty. At the time of this study, eight Annex 2 states have still not ratified the treaty, including North Korea. To detect, locate, and characterize nuclear explosions, an International Monitoring System (IMS) was established by the CTBTO as part of a verification regime. The IMS features four different approaches for the monitoring of potential nuclear explosions. Three methods (seismology, infrasound, and hydroacoustics) are attributed to waveform technologies and have the purpose of detecting, localizing, and identifying suspicious events with an explosive source mechanism. The fourth approach features the monitoring of particulate radionuclides and noble gases in the atmosphere and potentially provides unambiguous evidence of the nuclear character of an explosion. To assess the connection between radionuclide detections and potential source events, atmospheric transport modeling (ATM) is applied. Grouping data from different IMS technologies that is attributed to the same critical event is referred to as data fusion. For further information on the CTBT, the IMS, and the German National Data Center, the reader is referred to Pilger et al. (2017) or CTBTO (2018).

Since the first known nuclear test, carried out by the Unites States in 1945, more than 2000 confirmed nuclear explosion tests were conducted by China, France, Great Britain, India, North Korea, Pakistan, Russia, and the United States. It is generally accepted that all six past explosions in North Korea between 2006 and 2017 have been nuclear underground tests. This study therefore refers to the events in North Korea as nuclear tests, even though definitive proof may be missing. North Korea is the only country breaking the de facto moratorium on nuclear tests. For detailed information on the first five North Korean tests, see for example Hartmann et al. (2017).

All past North Korean nuclear tests were conducted at the Punggye-ri test site in the vicinity of Mt. Mantap in the northeastern part of the country. Following the geologic description of the test site area provided by Coblentz and Pabian (2015) and Pabian and Coblentz (2017), Mt. Mantap is made up of two distinct geologic formations. The core of the mountain consists of igneous basement rock of either diorite or granite, while the top is capped by a thin layer of basaltic lava flows. The basement crystalline rock and the top layer of Mt. Mantap are separated by a nearly horizontal sequence of volcanic deposits with a thickness of around 200 m. The origin of these deposits is suspected to be volcanic ash from Mt. Paektu, a volcano located around 100 km northeast. Due to their loose consolidation, the volcanic deposits are softer then the basement rock or the basaltic layer and are therefore more susceptible to erosion. Furthermore, the volcanic layer is more exposed to erosional scars and landslides due to its steeper slope. The erosion of the volcanic layer can cause the overlaying basalt cap to break off at the scarp, which is visible in larger piles below the volcanic layer where the slope of the mountain decreases slightly.

The official channels of North Korea announced the 2017 event as a successful test of a fusion bomb (Korean Central News Agency, 2017). This would be a major step in the nuclear program of North Korea. From a scientific point of view, therefore, the depth of the event, its strength in terms of radiated high- and low-frequency seismic energy, the contribution of possible faulting or slope instability processes, the near-surface damage in the test area, and the proof of whether fission products are detected as atmospheric tracers are key questions to be answered. These questions are approached by an integrated study based on different seismological (Sect. 2), infrasound (Sect. 3), remote-sensing (Sect. 4), radionuclide monitoring (Sect. 5), and modeling techniques which complement each other. As there is no easy concept to verify and characterize nuclear explosions, especially in a country where direct observations are difficult to assess, this study demonstrates and emphasizes the strength of an integrated multi-technology approach.

In this context, new methodical approaches are introduced to improve the depth estimate from teleseismic observations and to quantify uncertainties in the non-isotropic source component. The seismological study retrieves an independent absolute location based on a combination of seismological and remote-sensing data. The relative location between the six North Korean nuclear tests is obtained by means of waveform cross-correlation time lag data. The event depth is estimated for the first time by a joint inversion of source time function (STF) and waveform composition of direct and depth phases observed at small aperture, high-frequency arrays at teleseimsic distances. A full waveform moment tensor inversion (MTI) is applied and compared to results of previous explosions in North Korea. The source time overshoot and peak amplitude are compared to traditional body wave magnitude (m_b) estimates. Additionally, the effect of topography on peak amplitudes is estimated by 2-D-waveform modeling to assess the possible range of explosive yield. Infrasound observations and modeling are used to understand the earth–atmosphere coupling and the propagation of infrasound from the North Korean test site. The analysis of satellite-based remote-sensing data is important to improve the absolute location of shallow sources as well as to quantify the secondary mass movement effects at the surface. Radionuclide monitoring demonstrates the importance of ATM to avoid overinterpretation of variations in ¹³³Xe concentrations. While radionuclide monitoring provides the only direct evidence of nuclear explosions, it is demonstrated by careful modeling how difficult it is to interpret such data and that early claims of causal anomalies have possibly been over-interpreted.

2 Seismological investigations

2.1 Epicenter location

The localization in the Reviewed Event Bulletin (REB) of the CTBTO International Data Center (IDC) uses 125 seismic stations of the IMS and results in an epicenter of 41.321^∘ N and 129.035^∘ E with an error ellipse area of 110 km². This relatively high error is generated by the fact that only eight IMS stations are located at distances of 400 up to 2100 km (see Fig. 1a).

Figure 1 (a) Seismic stations of the International Monitoring System and other earthquake monitoring networks within a radius of 2100 km from the North Korean Punggye-ri nuclear test site. Stations that are used in the moment tensor analysis in Sect. 2.3 are marked with a small yellow dot. (b) Zoom into the Punggye-ri test site area. Numbered circles indicate the absolute locations of the six North Korean nuclear tests. The relative location error for each test is less than 100 m. The dashed line marks the profile used in Sect. 2.5 to study topographic influence of the Mt. Mantap on seismic energy generation and radiation from the test.

However, the incorporation of 25 additional seismic stations at regional distances cannot significantly improve the absolute location estimate, as the closest station MDJ is still 372 km away from the test site. Due to these large source–receiver distances, a further improvement of the absolute location accuracy is limited. Nevertheless, relative location methods can be applied for a high-precision localization of the events (Gibbons et al., 2017; Zhang and Wen, 2013; Zhao et al., 2014, 2016). In this study a relative location procedure based on the cross correlation of seismograms from 33 regional seismic stations for the six North Korean tests is applied. Seismograms for the past North Korean events significantly vary between the stations, due to propagation paths from source to the surrounding stations at different distances and azimuths. However, the correlation of the seismograms of the six events for each individual station shows a good coherence between the corresponding signals, with the Pn phase being the most pronounced arrival in the waveforms. These signals are correlated for each individual station, which has recorded at least two North Korean nuclear tests. The cross correlation between each pair of events is performed using normalization to 1.0 for the autocorrelation of each signal at zero lag. Maximum correlation values of 0.7 to 0.99 are obtained for stations of up to 1100 km distance for the four tests in the years 2009, 2013, and 2016. Due to significant differences in explosive yield, and therefore different source time durations, the 2006 and 2017 tests are slightly less correlated with the other tests. See Fig. 2 for the complete results of the cross-correlation analysis.

Figure 2 Maximum correlation coefficients for the coherence of the Pn phase for all available pairs of nuclear tests at 18 stations at regional distances of up to 2100 km. Due to their similar size in explosive yield, the four tests in the years 2009, 2013, and 2016 have a similar source mechanism and consequently show high maximum correlation coefficients between 0.70 and 0.99. The seismic signals for the 2006 and 2017 test are in general less correlated with the other tests due to the different source durations suspected for these two events. The 2006 and 2017 tests are the smallest and largest explosion, respectively.

The estimation of time lags turned out to be reliable, when the maximum cross-correlation values exceeded a threshold of 0.4. With this constraint, a subset of 165 event pairs for 19 stations was selected, for which precise travel-time differences with an accuracy on the order of the sampling rate used of 0.025 s could be determined. Pn-phase onset times of all tests are aligned to the onset at the closest station MDJ and fixed as relative start time for the estimation of the correlation time lags at the other stations. A double-difference method (Waldhauser and Ellsworth, 2000) is applied to cross-correlation time lags and results suggest that the last five events are located within a radius of 400 m, while the 2006 test is located around 2 km further to the east (see Fig. 1b). The relative locations of the six tests can be associated with absolute coordinates as soon as the geographical coordinates of one of the tests are known. The geographic, absolute location of the January 2016 nuclear explosion is fixed by means of radar interferometry data to the location of the maximum surface deformation observed after the test (Hartmann et al., 2017; Wei, 2017). The absolute epicenter location of the 2017 test from relative location procedures is consequently determined to be 41.3007^∘ N, 129.0728^∘ E (Fig. 1b).

2.2 Estimation of hypocenter depth and seismic moment

The source depth is needed to estimate the explosion process and strength of the seismic source. However, constraining the depth of a shallow source is difficult from regional and teleseismic data without a close station within a focal depth's distance from the source. The modeling of the wavelet consisting of a surface-reflected P phase, the so-called depth phase, and STF can potentially help in such a case. Such a wavelet modeling approach resembles the inverted and time-shifted STF and needs high-frequency waveforms above 1 Hz to resolve the onset of the depth phase, and it may be difficult to perform this from single station recordings if the signal-to-noise ratio (SNR) is poor at teleseimsic distances. An approach uncommon for nuclear test studies is used, which was established for the analysis of induced seismicity, where waveform beams are calculated at several small-aperture, short-period arrays to enhance the SNR (Fig. 3a).

Figure 3 (a) Setup of the depth analysis with the 2017 North Korean nuclear test marked as a red star and the four small-aperture, short-period arrays marked as bluish triangles. The epicentral distances to arrays IMAR, BCA, YKA, and TXAR are 48^∘, 54^∘, 65^∘, and 94^∘, respectively. (b) Residuals between observed and synthetic waveforms (L₂ norm) as a function of depth of the explosion. (c) Least squares source time function associated with the best-fit regularized solution. Solid and dashed lines are the normalized moment rate and moment function, respectively. (d) Observed (black) and synthetic (red) waveforms for the best-fit regularized solution at 600 m depth. The gray and green dashed lines indicate the theoretical arrival times of the P and pP phases, respectively. The traces are bandpass-filtered between 0.5 and 2.5 Hz.

The beam waveforms represent a superposition of direct and reflected waves. The depth of the explosion is estimated by comparing observed and synthetic beam waveforms (Fig. 3d), which are calculated for the best moment tensor solution described in Sect. 2.3, a common STF and varying depths. The source time is represented as a composition of multiple basis functions with unknown weighting coefficients, which are estimated in a least squares inversion with a smoothing constraint (Fig. 4). If u(t) is the seismogram (beam) at the array with coordinates r_r and from a seismic source with coordinates r_s, the P-wave train (i.e., P plus pP plus sP phase) at time t is represented by

$$\begin{array}{}\text{(1)}& u\left(t\right)=\left\{M_{jk}{G}_{j,k}({\mathbf{r}}_s,{\mathbf{r}}_r,t)\right\}\cdot m\left(t\right),\end{array}$$

where ⋅ is a time convolution and where the summation convention is applied. G_j,k are the spatial derivatives of the Green's function, where the comma before index _,k indicates a spatial derivative with respect to x_k. It is assumed that all moment tensor components M_jk have the same time dependency, which is described as normalized moment function m(t) with $m(t\to \mathrm{\infty })=\mathrm{1}$. The displacement in the far field is controlled by the time derivative of m(t), which is declared as moment rate function $\dot{m}\left(t\right)$. This explains that the far-field body wave pulses from earthquakes are single-sided pulses (Dahm and Krüger, 2014). Since only the far-field wave field is considered, the Green's functions in Eq. (1) are replaced by far-field Green's functions G^(ff), and m(t) is replaced by $\dot{m}\left(t\right)$. The inversion is set up for $\dot{m}\left(t\right)$ using a set of N time-shifted triangular basis functions h_l(t) (Fig. 4):

$$\begin{array}{}\text{(2)}& \dot{m}\left(t\right)=\sum _{l=\mathrm{1}}^N{\mathit{\beta }}_l{h}_l\left(t\right),\end{array}$$

where β_l is a weighting factor and each h_l(t) satisfies $\int h_l\left(t\right)\mathrm{d}t=\mathrm{1}$.

Figure 4 (a) Schematic sketch of five triangular basis functions h_l(t) with an equal weight of 1 and with a duration of 0.4 s each. The focal sphere symbols indicate that the orientation of the source mechanism is assumed to be constant during the rupture process. (b) Simulation of a complex moment rate function (blue line) by applying different weights to h_l and linear superposition.

Hence, Eq. (1) is written as

$$\begin{array}{ll}{\displaystyle }u\left(t\right)& {\displaystyle }=\left\{M_{jk}{G}_{j,k}^{\left(\mathrm{ff}\right)}\left(t\right)\right\}\cdot \left[\begin{array}{cccc}{h}_{\mathrm{1}}\left(t\right)& h_{\mathrm{2}}\left(t\right)& \mathrm{\dots }& h_N\left(t\right)\end{array}\right]\\ \text{(3)}& {\displaystyle }& {\displaystyle }\cdot \left[\begin{array}{c}{\mathit{\beta }}_{\mathrm{1}}\\ {\mathit{\beta }}_{\mathrm{2}}\\ \mathrm{⋮}\\ {\mathit{\beta }}_N\end{array}\right].\end{array}$$

The convolution in Eq. (3) can be written in discrete form, leading to an overdetermined matrix system for unknown weighting factors to be solved in a least squares sense using the L₂ norm (Dahm and Krüger, 2014). If more than one array is available, the equations are added to the coefficient matrix to realize a joint inversion. Additionally, the STF inversion is stabilized by using a regularization with a roughness matrix. The STF is inverted for every single trial depth. The length of the STF is constrained to 1.5 s, where 0.5 s is considered before the arrival of the P wave. It should be noted that, with a shorter length of STF, the details of the STF with a complex structure are difficult to resolve. On the other hand, with longer length of the STF, i.e., a larger number of basis functions, the fitted curve oscillates wildly and gives a very poor representation of the model; this is known as over-fitting. The length of the STF mentioned before was chosen to reach a compromise, after examining the residual misfit in total signal variance computed for different lengths of source time function. In order to compensate for travel-time residuals occurring along the mantle ray, zero-lag Green's functions are correlated with P-wave beams at every array and the travel-time corrections are considered before the inversion. The source depth itself is estimated using a grid search approach, where the depth of the point source Green's functions is sampled from 0 to 1500 m in 100 m steps. For every depth solution the residuals are estimated and stored for final evaluation.

Green's functions are calculated using the code QSEIS (Wang, 1999), where source- and station-specific crustal models can be considered. AK135 (Kennett et al., 1995) is used for the mantle, while different crustal profiles for source and station regions are taken from CRUST 2.0 (Bassin et al., 2000). Intrinsic attenuation for P waves is set to 5000, since otherwise high frequencies are damped out at teleseismic distances. The sampling frequency is 20 Hz. The grid search depth phase modeling has been applied previously to different cases of induced seismicity (Dahm et al., 2007) and the 2016 nuclear explosion (Cesca et al., 2017), but the simultaneous STF inversion is implemented for the first time in this study.

The full waveforms of beams, their peak heights, and peak distances are well fitted in a joint inversion of four arrays. The best fits are obtained at depths between 400 and 800 m (Fig. 3b). However, it should be noted that the depth value obtained from the inversion procedure might be overestimated due to the use of an unperturbed velocity model. Perturbations in the velocity model, for example through potential fracturing above the explosion, which would substantially reduce P-wave velocity, were neglected. If the source of the synthetic beams is placed deeper or shallower, the duration and appearance of the P , pP, and sP pulses changes and residuals increase. Of interest are the retrieved moment rate and moment function, which show a clear double pulse and overshoot, respectively (Fig. 3c). Such an overshoot is not expected for the rupture process of tectonic earthquakes, where the moment rate and moment functions of the P wave are single-sided pulses and monotonously increasing functions, respectively. For nuclear explosions, however, this feature is commonly observed. This can be explained by at least a partial collapse of the explosion cavity immediately after the explosion. In this case, the final moment is only 23 % of the peak moment in Fig. 3c. The inversion of long-period waves led to a seismic moment M₀ of 2.33×10¹⁷ Nm, representing the time after the overshoot. The peak seismic moment is thus estimated in the range of $M_{\mathrm{0},\mathrm{peak}}\approx \mathrm{1.02}\times {\mathrm{10}}^{\mathrm{18}}$ Nm. The associated low- and high-frequency moment magnitudes (M_w) are M_w=5.55 and $M_{\mathrm{w}}^{\left(\mathrm{peak}\right)}=\mathrm{5.97}$, respectively, and can explain the large difference between the long-period M_w and the high-frequency m_b estimates which emerged during the magnitude estimation described in Sect. 2.4.

2.3 Moment tensor inversion of the test and the main aftershock

Seismic signals of the 2017 North Korean test show a great similarity to those generated by the previous four nuclear tests conducted in 2009, 2013, and 2016. Especially long-period waveforms with periods above 10 s recorded at regional distances show a high waveform similarity and a substantially increased amplitude, revealing a very similar radiation pattern of all events and a larger moment release for the 2017 nuclear test. In recent years, MTIs have been performed for nuclear tests in North Korea, including those carried out in 2009, 2013, and 2016 (Barth, 2014; Cesca et al., 2017; Ford et al., 2010; Hartmann et al., 2017; Vavryc̆uk and Kim, 2014). All MTIs were performed by fitting low-frequency full waveform seismic data at regional distances, either in the time or in the frequency domain. The majority of these solutions revealed a significant positive isotropic component. For the 2013 nuclear test an extraordinarily high double couple (DC) component was found, indicating differences in containment or near-source damaging effects (Barth, 2014; Vavryc̆uk and Kim, 2014). For the 2017 test the seismic moment tensor is inverted by fitting the low-frequency amplitude spectra (epicentral distances up to 1200 km) and full displacement waveforms (epicentral distances up to 600 km) in the frequency range from 0.02 to 0.04 Hz following the approach by Cesca et al. (2013, 2017). For the moment tensor optimization, the Grond algorithm (Heimann et al., 2018) is used, which searches the parameter space by random search first and later by an increasingly denser search around the best solutions. The algorithm finally resolves the full moment tensor components, the centroid location, and centroid time and also provides information on model uncertainties and source parameters trade-offs. Full waveforms and their spectra are calculated assuming a layered crustal model proposed by Ford et al. (2010). Examples of full waveform displacement and amplitude spectra fits are shown in Fig. 5.

Figure 5 Full waveform displacement and amplitude spectra fits (gray and red lines denote observations and synthetics, respectively) are shown for the nuclear explosion for the vertical (a), transversal (b), and radial components (c). The top left side of each plot reports the trace component (network, station, location, and component codes), epicentral distance, and azimuth. Waveform plots report the starting time (relative to the origin time) and duration (in minutes and seconds) of the time window used (the gray background represents the applied taper), while spectra plots report the frequency range in Hz. For geographical locations of the stations used, see Fig. 1a.

The moment tensor solution for the 2017 test shows a dominant positive isotropic part of 60 %, 16 % positive compensated linear vector dipole (CLVD), and 24 % of DC. The best solution, found at a depth of around 2 km, has a scalar moment of 2.33×10¹⁷ Nm, equivalent to a M_w of 5.55. The centroid depth is poorly resolved and good waveform fits are found for very shallow sources down to 2.5 km. This result confirms the very shallow depth of 400–800 m accurately resolved by the analysis of depth phases in Sect. 2.2. The moment tensor optimization resolves a broad ensemble of well-fitting moment tensor solutions. They are plotted using a source type diagram in Fig. 6 and are characterized by strong deviations from a pure DC model (mapping in the center of the diagram), a typical model for tectonic earthquakes. Instead, acceptable solutions are either dominated by a positive isotropic component (positive values of the y axis in the plot), a negative vertical CLVD (positive values of the x axis), or a combination of both. The variability among these type of solutions has been attributed to trade-offs among moment tensor components for very shallow sources (Cesca and Heimann, 2017; Cesca et al., 2017). This range of moment tensor solutions includes alternative moment tensor configurations as proposed for previous nuclear explosions in North Korea (Barth, 2014).

Figure 6 Best moment tensor decomposition (top) and source type plots (bottom) for the 2017 explosion. Top: the moment tensor decompositions reports the focal sphere for the full moment tensor, and its isotropic and deviatoric components. The deviatoric component is further decomposed into CLVD and DC terms. The relative size of focal spheres scales with their relative magnitude. Bottom: the source type diagram after Hudson et al. (1989) is used to judge the moment tensor decomposition, illustrating the ensemble of best-fitting solutions. Figure shows both solutions based on all available data, as well as using data subsamples upon bootstrapping. The position of each solution (circles, focal sphere for selected solutions) illustrates the moment tensor (MT) decomposition, while the focal spheres (shown for selected solutions only) illustrate the geometry of the DC term. The decomposition of the best solution is denoted by a blue square.

In the aftermath of the nuclear test, a seismic event took place in its direct vicinity around 8 min later. The estimated scalar moment and M_w are 1.88×10¹⁶ Nm and 4.81, respectively, and thus the scalar moment is only ∼8 % of the one of the nuclear test. The waveforms of the aftershock still show a sufficient SNR to perform an MTI. Waveform displacement and amplitude spectra fits of the MTI are shown in Fig. 7.

Figure 7 Same as Fig. 5 but for the aftershock around 8 min after the test.

The inversion is based on spectral and waveform fits at stations located within 1000 km epicentral distance, filtered between 0.02 and 0.04 Hz. The calculated depth of the aftershock is similar to the depth of the explosion. The best MTI solution suggests a dominant implosive source component (65 %), negative CLVD (29 %), and a DC component of 6 %. This source mechanism is compatible with a shallow collapse source or might hint at some kind of break-in process. An ensemble of well-fitting moment tensor solutions for this inversion is shown in Fig. 8.

Figure 8 Same as Fig. 6 but for the aftershock around 8 min after the underground test.

2.4 Magnitude and yield estimation

The determination of m_b is the prerequisite for a reliable estimation of the yield of a nuclear explosion. In order to obtain comparable results for the past six North Korean events, m_b is calculated from waveform peak values measured at always the same 15 IMS seismic stations that recorded all six events. These stations (MKAR, INK, FINES, WRA, YKA, ASAR, AKASG, HFS, NOA, VRAC, GERES, STKA, DAVOX, NVAR, PDAR) cover the entire azimuthal range with a maximum gap of 121^∘ and are located at distances between 3700 and 9000 km from the Punggye-ri test site. Exemplarily, the seismic signals at the German primary seismic station GERES are shown in Fig. 9. The high similarity between the six nuclear explosions stands out due to the short and simple P-wave train. The recordings differ only in their amplitudes. The resultant single station magnitude for GERES was used together with the corresponding magnitudes of the other 14 stations involved to obtain reliable and comparable average magnitudes.

Figure 9 Beam seismograms calculated from the vertical components of the GERES array showing the P wave of the six North Korean nuclear explosions. The seismograms are normalized to the individual maximum amplitudes, for which the station magnitude is estimated.

For the 2017 test a m_b of 6.2 is calculated, compared to m_b values of 4.1 to 5.3 for the five preceding tests. In general it is not possible to state one single relation between magnitude and yield, as this relation depends on various factors, such as the geological setting at the source site, the efficiency of wave propagation from source to receiver, the depth of the explosion, and the coupling of the source to the underground. A number of empirical formulas of the type $m_{\mathrm{b}}=A+B\log \left(Y\right)$, with Y being seismic yield in kT TNT equivalent and A and B being constants depending on the aforementioned factors, have been developed to relate yield and m_b. These relations have been successfully used for yield calibration, for example at the Nevada test site (Murphy, 1981), in Kazakhstan (Ringdal et al., 1992), or in Nova Zemlya (Bowers et al., 2001). As no particular magnitude–yield relation has been approved so far for the North Korean test area, the latter relation by Bowers et al. (2001) is used in this study. The North Korean nuclear explosions were mainly conducted beneath the Mt. Mantap complex, consisting of stratified volcanics covering a basement of granite and diorite. In general, wet hard rock can be taken into account for these geological conditions, when the yield is estimated on the basis of the seismic magnitude. Under these assumptions a yield of around 400 kT TNT equivalent is estimated for the September 2017 test, clearly illustrating a steady increase in explosive strength of the nuclear tests (see Fig. 10).

Figure 10 Magnitude–yield relation curves for different geological settings. Numbered circles indicate the six North Korean nuclear tests. Gray background shading represents lower and upper boundary literature values for the different environments.

2.5 Influence of the Mt. Mantap topography

Two-dimensional synthetic wave-field simulations are used to study the influence of topography of the Mt. Mantap complex on the propagation of P- and S-wave energy emitted by a point-like explosion source. A Chebyshev pseudo-spectral method (Tessmer et al., 1992) is used for modeling the elastic wave fields on a Cartesian 820×512 grid for an 8 km × 5 km box in x and z direction, respectively. No variation in medium properties is assumed in the y direction. Figure 11 shows the model configurations as well as the wave propagation for an explosion source at 0.6 km depth below the surface (see Sect. 2.2).

Figure 11 Snapshots obtained by and velocity model used for numerical simulations. P-wave velocities are indicated by shading according to the gray scale map (for S-wave velocities we assume a v_P∕v_S ratio of 1.73), as well as logarithmically scaled divergence (P-wave energy) and curl (SV-wave energy) by bluish and greenish colors, respectively. The snapshots in (a) show the P- and SV-wave separated propagation for an explosion source without topography, the middle panels with topography across the North Korean test site for a west–east profile along 41.3^∘ N with Mt. Mantap in the middle. Panels in (c) are similar to those in (b) but with a more realistic geological setting, where erosion is accounted for. The time stamp of each snapshot is displayed in the upper left corner. The point source is beneath the center of Mt. Mantap at a depth of 0.6 km below the surface. The panel below the snapshots compares the summed average divergence and curl computed every 0.02 s in a rectangular box below the source for the case without (yellow) and the cases with (green and blue) topography. The downward-propagating wave field shows strongly increased amplitudes for pP and pS phases in both cases of topography.

The source time function is taken from Glasstone and Dolan (1977), where the source duration is 0.15 s. The wave field V is separated in divergence ($|\mathrm{\nabla }\cdot V|$) and curl ($|\mathrm{\nabla }\times V|$) reflecting the P and SV energy, respectively. Moreover, the wave propagation is considered for one model without and two models with topography, whereas the topography is considered for a west–east profile along 41.3^∘ N crossing the Mt. Mantap massif (see Fig. 1b). The velocity model is chosen such that it reflects the petrology of the local geology as described by Coblentz and Pabian (2015). However, strong lateral variations as thrust or normal faulting in the uppermost layers are not taken into account due to reasons of numerical stability. Instead lateral variations, which are based on an exponential perturbation with a correlation length of 100 m, are considered with differences in density, P-, and S-wave velocities of ±3 %. Merely the last model accounts for erosion of the uppermost basaltic and volcanic layers as well as an almost vertical dyke reflecting the transition between the granodiorites and the granite from west to east. Nevertheless, computing the P divergence and SV curl in a rectangular box shows large differences in the downward-propagating wave fields, where strong increases for both pP- and pS-wave amplitudes can be observed. The position and the lateral extent of this box covers a slowness range from −9 to +9 s ^∘−1, which are typical slowness values for teleseismic P waves. For a source under Mt. Mantap, surface reflections are increased in their amplitudes by a factor of ∼1.5 due to topography.

Wave energy is focused due to the shape of the mountain, especially from the easterly flank with a change in altitude of approximately 600 m (see Fig. 1b). Topography is the crucial factor; the difference between the two simulations with and without a realistic geologic setting (see Fig. 11b and c) is negligible. Further simulations for sources shallower and deeper than the reference value of a depth of 0.6 km below the surface provided similar pictures, whereas for sources shallower than 0.5 km and deeper than 1.7 km, respectively, energy focusing is less pronounced with amplifications in amplitude of 1.3 and less. Although only a two-dimensional model is considered, such simulations are also valid for the three-dimensional space because Mt. Mantap exhibits a high symmetry. In general, this numerical modeling gives indications of (1) the generation of clear infrasonic signals (see Sect. 3) because topography (slope of Mt. Mantap) increases the amount of radiated acoustic energy in the subhorizontal direction, which is the usual propagation path for infrasound (compare Fig. 12); (2) a reduced value of the isotropic part obtained by MTI because a larger amount of S-wave energy is generated in the source region due to topography effects above the source (see Sect. 2.3); (3) an overestimation of the yield of the explosion because P and pP phases are considered as a single onset in teleseismic distances and in the investigated frequency range for shallow sources, leading to an increase in maximum P-wave amplitude of approximately 1.3 and hence in m_b of around 0.1; (4) no effect on differential travel times between P and pP phases due to topography and therefore no effect on hypocenter depth estimation using depth phase modeling (see Sect. 2.2).

3 Infrasound observations

For a shallow event in the given magnitude range, the observation of infrasound signals at distances of up to hundreds of kilometers is generally expected (Mutschlecner and Whitaker, 2005). Infrasound is subaudible sound below the human hearing threshold of 16 Hz and is the dedicated monitoring technology of the CTBTO for remotely detecting atmospheric explosions. Furthermore, infrasound can be used in civil and scientific applications such as the monitoring of volcanoes, meteoroids, ocean swell, and various other natural and anthropogenic activities (see Le Pichon et al., 2010 for a comprehensive overview). Infrasound signals are also generated by strong underground or underwater sources when the generated waves couple into the atmosphere. The North Korean tests between 2006 and 2016 were investigated and registered by various infrasound arrays (Assink et al., 2016; Che et al., 2014; Koch and Pilger, 2018; Park et al., 2018). Detections of the 2013 test by the IMS stations I30JP in Japan and I45RU in Russia were part of the REB and thus confirmed the underground nuclear test as a source of infrasound waves. For the 2017 test, the station I45RU, located at a distance of 401 km from the test site, detected infrasound signatures from the nuclear test (Assink et al., 2018). Clear waveforms undoubtedly associated with the nuclear test were recorded and were again part of the REB. Figure 12a highlights the waveform beam of the Russian four-element infrasound array I45RU (denoted as the colocated seismic station USRK in Fig. 1a).

Figure 12 (a) Waveform beam of the four-element infrasound array I45RU with corresponding celerity values. (b) Progressive Multi-Channel Correlation analysis in the same time frame for signal frequencies between 0.2 and 4 Hz. Information on back-azimuth direction is color coded. True back azimuth from test site to station is 218^∘. (c) Propagation modeling using ray tracing with gray-shaded travel-time information for eigen rays from the test site to the station and color-coded attenuation information from parabolic-equation modeling. Is₂ and It denote the corresponding stratospheric and thermospheric arrivals after 1450 and 1520 s, respectively. Is_1f shows the early signal detection and theoretical ray path for an infrasonic forerunner.

Coherent waveform signals can be observed and associated with infrasound propagation from the nuclear test site arriving at I45RU after travel times of 1450 and 1520 s, corresponding to celerities of 290 and 270 m s⁻¹, respectively. Preceding waveform activity at around 1350 s (315 m s⁻¹) can either be associated with infrasound radiation from an extended (near-)epicentral source region or the fast acoustic phase of an infrasonic forerunner (Evers and Haak, 2007). These values are in the pure acoustic range and indicate stratospheric ducting (Is₂ with two reflections and Is_1f with one elongated reflection, turning heights of 40 to 55 km) and thermospheric ducting (It, turning heights above 100 km) (Drob et al., 2003). It is the first time that an underground nuclear test is strong enough to generate distinct and far-reaching pressure signatures that were clearly recorded at a remote IMS infrasound station after propagation through the thermosphere above 100 km altitude. The corresponding Progressive Multi-Channel Correlation (Cansi, 1995) analysis of the signal content in terms of time variation, frequency content, and back-azimuth direction is illustrated in Fig. 12b. The PMCC method draws rectangular time–frequency windows wherever correlated signal energy in three or more array elements is present. Color-coded back azimuths between 210 and 225^∘ (southwest) point towards the Punggye-ri nuclear test site with a true back azimuth of 218^∘. Effects from cross winds (west to northwest, strongest in the stratosphere) lead to the fact that the observed back azimuths are some degrees lower, especially within the stratospheric duct. The PMCC signatures with broadband frequency content from 0.2 to 4 Hz can thus be clearly associated with the 2017 North Korean nuclear test.

Figure 12c models the propagation from the test site epicenter to the station I45RU using a parabolic equation and ray-tracing methods for the attenuation of signal amplitude and the connection of source and receiver by eigen rays. A combination of European Centre for Medium-Range Weather Forecasts (ECMWF) analysis data and climatological profiles is implemented for temperature and wind model backgrounds including gravity wave perturbations (see Koch and Pilger, 2018, for further details on the methods and the atmospheric profiles used). The eigen rays show that the stratospheric phases are about 170 and 70 s faster than the thermospheric phase (caused by a higher effective sound speed in the stratosphere and the longer wave path for the thermospheric phase). The signal attenuation, as obtained from parabolic-equation modeling, indicates that only a small portion of signal energy is ducted in the stratosphere caused by partial reflections from gravity wave variations in the stratospheric mean background. This leads to about a 4 to 5 dB higher attenuation in the stratospheric duct and thus stronger signal amplitudes in the thermospheric duct by a factor of 1.6 to 1.8, which corresponds to the observed waveforms. The presented infrasound signals clearly observed even at 401 km distance and with disadvantageous propagation conditions are indicative of a strong surface movement in the epicentral area (see Sect. 4) of the underground test generating large acoustic pressure peaks as argued in Sect. 2.5. Apart from the strong epicentral surface movement, infrasonic signatures were also identified from seismo-acoustic coupling and the assumed cavity collapse associated with the 8 min subsequent aftershock (see the Supplement to this paper for further information).

4 Remote-sensing studies

To visualize and characterize the surface imprint of the 2017 test, radar data from the ALOS-2 satellite (Suzuki et al., 2009) and multispectral optical data from the Pleiades satellite (Gleyzes et al., 2012) are investigated. Data from the ALOS-2 satellite in the L band (1257.5 MHz frequency, 3 m resolution) are analyzed by interferometric synthetic aperture radar (InSAR) to investigate surface deformations (subsidence and uplift) correlated with the 2017 test and its aftershocks. The method of change detection analysis based on the Normalized Differential Vegetation Index (NDVI) is applied. Spaceborne InSAR has been used to detect surface displacements of the Earth's surface for more than 20 years (Massonnet et al., 1993) and has reached maturity in recent years (Adam et al., 2009), nowadays allowing for the detection of surface displacements with a precision on the order of a few millimeters. Repeat pass interferometry is based on the registered interferometric phase per ground cell (pixel) and is related to the distance differences between the scatterer and the synthetic aperture radar (SAR) sensor between two acquisitions separated in time. The interferometric phase measured as moduli 2π is thus ambiguous. To solve this ambiguity, frequency estimations using unwrapping techniques are calculated. The Earth's curvature, topography, and atmospheric effects influence the interferometric phase in addition to surface displacements. These phase contributions have to be estimated and subtracted to correctly estimate the surface displacement (Kampes, 2006). In addition, the backscattered electromagnetic signal between the acquisitions times needs to be correlated or coherent. Decorrelation can be caused by (1) temporal (large separation between the acquisitions in time), (2) geometrical (large perpendicular baseline) and (3) scattering effects (Gatelli et al., 1994).

Differential InSAR (DInSAR) processing was adjusted to the test site target area regarding the abovementioned constraints. As the investigated area is located in a mountainous region, DInSAR processing is challenging due to factors such as limited visibility for SAR acquisition geometry, tropospheric phase contributions, and/or snow cover and vegetation causing temporal phase de-correlations. DInSAR processing encompasses the following steps: co-registration, filtering, unwrapping, and the discussion of coherence. Coherence serves as a measure of the quality in DInSAR studies, with coherence values ranging from 0.0 to 1.0, with 1.0 suggesting a perfect interferogram with no noise. For the analysis of the surface displacement due to the 2017 test, data from 29 August and from 12 September 2017 are used, data coherence between these two times has a mean value of 0.6 for the investigated area. For areas that show surface deformation, the coherence values range from 0.2 to 0.8.

Figure 13 shows the surface deformation at the test site during the 2017 test restricted to pixels with coherence values greater than 0.25.

Figure 13 Surface displacements from DInSAR analysis for the 2017 North Korean nuclear test in September 2017. Coherence threshold is set to 0.25. Data from the days 29 August and from 12 September 2017 are used. Arrows indicate slope direction and angle. Numbered circles indicate the epicenters of the tests from 2006 to 2017. Black solid box delineates an area of around 3×4 km² in which the surface displacement is between −10 and 10 cm.

In the central part of the Mt. Mantap massive, coherence values dropped significantly below 0.25 due to very strong ground deformations (Wang et al., 2018). The recognition of the displacements is only possible in the line of sight in the direction towards or away from the sensor. Due to the incidence and look angle of the sensor measured with respect to the ground and north direction, respectively, it is not possible to fully resolve the amount of the vertical surface displacement caused by the 2017 nuclear test. During data acquisition on 29 August and 12 September 2017, the sensor had an incidence angle of 43^∘ and a look angle facing east–northeast. These angles combined with the topography of Mt. Mantap result in the maximum measurable percentage of vertical surface displacement. For the western flank of Mt. Mantap pointing towards the sensor with a slope of around 27^∘, the percentage of resolvable surface displacement is only around 30 %; for the eastern flank pointing away from the sensor with a slope of around 10^∘, the value of maximum measurable surface displacement increases to around 80 %. The resulting map of resolvable displacements clearly shows an area of subsidence of up to 10 cm around 3 km north of the main support area (compare Fig. 1b) and an area with clear uplift of up to 10 cm west of the Mt. Mantap peak. DInSAR processing of C-Band Sentinel data and Terra SAR-X data for the 2017 test did not reveal acceptable quality measures regarding coherence in the central part of Mt. Mantap, whereas for the tests in 2016 Sentinel data did show excellent results (Hartmann et al., 2017; Wei, 2017). A study by Wang et al. (2018) shows additional results based on amplitude images from Terra SAR-X data and calculates 0.5 m subsidence in the center and 3.5 m of horizontal displacement. This strong spatial and temporal coincidence of the displacement areas to the nuclear test site suggests a high correlation with the underground test. To validate the displacement maps of the DInSAR analysis of the nuclear test, Pleiades data sets from 26 August and from 8 September 2017 were processed to reveal surface characteristics related to the underground test (Fig. 14). This relation is concluded by the strong temporal and spatial connection of the landslides to the test and aftershock origin time and epicenter locations.

Figure 14 Change detection analysis via NDVI (Normalized Difference Vegetation Index) from Pleiades data showing potential landslides (in purple patches) triggered by the 2017 nuclear test. Numbered circles indicate the epicenters of the tests from 2006 to 2017.

Change detection analysis show numerous landslides activated during the test and aftershocks. As a result of the processing of the ALOS-2 data, an area of around 3×4 km² can be delineated, where surface movement rates range between −10 and 10 cm. For comparison, the area of the January 2016 North Korean nuclear test, which showed comparable surface movement rates, was in the range of 1.5×1.5 km² (Hartmann et al., 2017). This large difference clearly underlines the enormous power of the 2017 nuclear test. The strong surface deformations at the surface observed for the 2017 test might lead to pathways and the subsequent release of radionuclides, especially gaseous radionulcides such as ¹³³Xe, from the test site. Potential measurements of radionuclides are discussed in the following Sect. 5.

5 Radionuclide monitoring and atmospheric transport modeling

Measurements of radioactive fission products released into the atmosphere even at great distances can deliver definitive proof of the nuclear character of an explosion. To assess the consistency between measurements and potential source locations, the ATM method is applied (Becker et al., 2007; Ross et al., 2017). For the ATM analysis of potential detections of radionuclides the Lagrangian particle dispersion model HYSPLIT (Stein et al., 2015) is used. The ATM is operated with meteorological data from the US National Centers for Environmental Prediction (NCEP). For the prediction of potentially affected radionuclide measurement stations, HYSPLIT is operated in forward mode with hypothetical releases. For the assessment of the potential source region of detected radioisotopes, the model runs in backward mode. Dispersion simulations of a potential immediate release after the explosion on 3 September 2017 show prevailing wind directions to the northeast (Fig. 15a).

Figure 15 Atmospheric transport modeling results. (a) Forward dispersion modeling after 36 and 96 h for an immediate release on 3 September 2017. (b) Number of overlapping backward plumes for seven ¹³³Xe detecting samples in northeastern South Korea from 7–10 September 2017 (fixed to 42 h transport time for each sample). (c) High-resolution forward dispersion modeling after 36 h for a delayed leakage on 4 October 2017. (d) Number of overlapping backward plumes for five ¹³³Xe detecting samples at RN58 Ussurysk during October 2017 (transport time fixed at maximum sensitivity to North Korean test site).

According to the forward simulations, the IMS station RN58 (Ussurysk, Russia) would have been affected, but it was not operational at that time. Station RN38 (Takasaki, Japan) would have been missed by the plume from the potential release from the test site which is typical of the wind patterns of the season. Stations further downwind, for example in North America, would have required a larger release for activity concentrations to exceed the detection threshold.

The South Korean Nuclear Safety and Security Commission (NSSC) issued a press release on 13 September 2017 with results of national efforts to detect radioactive Xe isotopes after the test. The statement “Considering the timing and location of sampling and air currents, the NSSC could conclude that the air currents from the North Punggye-ri area has been mixed to the sample and the Xe-133 detected from this air sampling is relevant to the 6th nuclear test of DPRK” (Nuclear Safety and Security Commission, 2017) was interpreted by the media as a source attribution. Backward simulations for the samples collected at the northeastern stationary measurement system at Geojin were performed using 500 000 model particles per sample and 0.5^∘ GDAS (Global Data Assimilation System) meteorological data. Additive overlay of the source regions of seven stationary 12 h samples with elevated ¹³³Xe from 7 to 10 September 2017 points to the region where the Yongbyon North Korean Scientific Research Center is located, as a likely coincident source (Fig. 15b), whereas the test site can be disregarded as a common source for the seven samples. This contradicts the suggestions in the media that the detected ¹³³Xe originated from the test explosion.

Continued forward simulations for hypothetical releases from the test site show that the station RN58 would have been repeatedly affected by the plumes during October. Figure 15c shows as an example the dispersion after 36 h of a hypothetical 10 h release on 4 October 2017 using HYSPLIT with 0.25^∘ GFS (Global Forecast System) data and 1 million particles released. In fact, the station RN58 went back to operation in between, and three peaks of elevated ¹³³Xe activity concentration were measured on the days indicated by routine forward simulations from the test site (5, 19, and 26 October). Those peaks consist of five samples with ¹³³Xe activity concentrations between 0.5 and 1 mBq m⁻³. Backward ATM results also indicate a high sensitivity of these five samples to being released from the test site. Figure 15d shows the overlap of the corresponding retro-plumes at a time of maximum sensitivity to the test site. A leakage on the order of 10¹¹ Bq¹³³Xe would have been sufficient to explain the detections at station RN58. Although the measurements and the ATM are consistent with delayed small releases from the test site in October, they cannot be unambiguously attributed to the North Korean test explosion.

6 Discussion

A combination of seismological analysis using a double-difference method and remote-sensing technologies results in an epicenter location of the 2017 test of 41.3007^∘ N, 129.0728^∘ E, clearly placing the event inside the North Korean nuclear test site, at which all previous tests have been conducted. The epicenter location is confirmed by results from Wang et al. (2018), who estimated an epicenter location of 41.300^∘ N, 129.078^∘ E ±50 m, which places the two epicenter solutions only ∼400 to 500 m apart. On 24 May 2018 North Korea invited international media to witness the demolition of the Punggye-ri test site. During a demonstration North Korean officials presented a map showing the relative locations of the six nuclear tests (CNN, 2018). Although the correctness of the information provided by the North Korean government cannot be verified, comparisons indicate that the assumed locations fit the results from the independent relative location method reasonably well. Two independent methods showed that the 2017 North Korea event was very shallow. The full waveform moment tensor inversion indicates a centroid depth within the range of about 2 km but has a relatively poor resolution because of the long wavelength of low-pass-filtered data. The depth phase modeling of the P waveform was performed at 2.5 Hz and found a centroid depth between 400 and 800 m. According to the crustal profile at the source side, the expected P-wave wavelength is ∼1000 m, which indicates that the source depth is still resolved. This is supposedly the best absolute centroid depth estimate possible from seismological data far from the source. Depth estimations in this study are comparable to results from Wang et al. (2018). The authors used numerical models to minimize the misfit between predicted and observed surface displacements and specified a depth of 450±100 m.

The explosive character of the 2017 event is confirmed by MTI analysis, which indicates a dominant isotropic source part of around 60 %. These results are supported by studies by Liu et al. (2018) and Han et al. (2017), who also revealed a dominant positive isotropic source part. Under the assumption of the magnitude–yield relation following Bowers et al. (2001), the estimated m_b of 6.2 (see Sect. 2.4) for the event results in a yield of 400  kT TNT equivalent. This value is assumed to be an upper limit, as the following factors might lead to a yield overestimation. (1) Topography at the test site (see Sect. 2.5) strongly influences the radiation of seismic energy from the source. This results in higher amplitudes of the pP phase. As the hypocenter location is very shallow, P and pP phases from the event arrive nearly simultaneously at distant receivers and superimpose at frequencies relevant for the estimation of m_b. This superimposition leads to an increase of approximately 1.3 in P-wave amplitude in the far field and consequently to an increase in m_b of around 0.1. (2) Furthermore, in the standard calculations of m_b (Bormann, 2012; Gutenberg, 1945a, b), distant-dependent correction curves, which are valid for a DC source, are applied. Gaebler and Ceranna (2017) state that the use of correction terms for a shallow explosion source can lead to an overestimation of m_b of around 0.2. Taking both these effects into account, the originally estimated value for m_b of 6.2 had to be considered as an upper limit estimate and has to be corrected by a factor of up to 0.3. This correction leads to a lower estimate for m_b of 5.9, similar to $M_{\mathrm{w}}^{\left(\mathrm{peak}\right)}$ that is estimated from the source time function, which translates to an estimated yield of around 160 kT TNT equivalent. Therefore, indications are strong that yield might be overestimated, not only for the 2017 explosion but also for all five previous North Korean nuclear tests.

Following the main event, around 8 min later, a large aftershock of local magnitude m_L 4.1 was detected. MTI analysis of this event shows a clear negative isotropic part which leads to the conclusion that this event is caused by at least a partial collapse of the cavity formed by the main explosion. These findings are again supported Liu et al. (2018) and Han et al. (2017), who also find a clear negative isotropic source part for the aftershock. According to Liu et al. (2018), the aftershock can qualitatively be interpreted as a rapid destruction of an explosion-generated cracked rock chimney due to cavity collapse. Following this strong aftershock, nine more earthquakes with local magnitudes m_L ranging from 2.5 up to 3.2 were detected by the Korean Meteorological Agency (Korean Meterological Agency, 2018) in the region (41.2 to 41.4^∘ N, 129.0 to 129.2^∘ E) until April 2018. Aftershocks of such a magnitude range and quantity were not observed from the five previous tests from 2006 to 2016 but are not unusual for nuclear explosions of that size (Boucher et al., 1969).

In the aftermath of the 2017 event measurable surface activity was detected in the region of the Punggye-ri test site. Remote-sensing results suggest a distinct pattern of uplift (up to 10 cm) in the western part of the Mt. Mantap massif and subsidence (up to −10 cm) in the eastern part. The geological interpretation by Coblentz and Pabian (2015) and Pabian and Coblentz (2017) shows a north–south trend of various geological units ranging from carbonates, diorites, granites, and gneises in contact with a younger basalt (layered volcanics) at Mt. Mantap. Furthermore, the structural pattern in the region suggests a north–south-trending system of faults up to the main tunnel entrances, with a part of the fault system branching out to the northwest. Considering this geological background the deformation pattern can be interpreted as a (partial) uplift of the basement rocks in the western part. Eventually this is structurally controlled by the north–south trending reverse faults together with the unconformable capping of stratified volcanic sequences. This interpretation is confirmed by multispectral optical data from the Pleiades satellite which shows that the layered tuffs are strongly affected by intensive landslide activities.

Overall, accounting for all these activities (surface deformations, the collapse event, aftershocks, landslides) it can be argued, that the internal structure beneath Mt. Mantap might be strongly stressed and shattered due to the last six tests and might be rendered useless for further test explosions (Wang et al., 2018).

Despite these strong activities and the enormous release of seismic energy from the event, no immediate release of radionuclides was measured by the IMS. This is a strong indication that no direct pathways from the cavity to the surface for the propagation of radionuclides were created during the test and the subsequent collapse-type aftershock. Only a few days after the 2017 event, South Korea announced national measurements of ¹³³Xe and proposed a connection with these observations to the test. However, backward ATM calculations performed in the context of this study disfavor the Punggye-ri test site as possible origin for ¹³³Xe but rather point to the Yongbyon nuclear complex as a most likely coincident source region instead in September. Nevertheless, in October 2017 three peaks of ¹³³Xe were detected at the IMS station RN58 in Ussurysk, Russia. Forward as well as backward ATM show directly connecting dispersion conditions and indicate the North Korean test site as a likely source for the ¹³³Xe detections. Nevertheless, other regional sources cannot be excluded, and it was not possible to definitely identify the presence of ^131mXe in the RN58 samples, which – in the correct isotropic ratio – would have increased the evidence of a nuclear explosion source (Kalinowski et al., 2010), as for example shown for the 2013 North Korean nuclear explosion by Ringbom et al. (2014). According to calculations of concentration and sensitivity, a release at the test site on the order of a few 10¹¹ Bq would have been sufficient to cause and explain the detections at station RN58. This is also quantitatively consistent as the required source term ranges well below the remaining fission inventory even after 6 to 10 half-lifes of ¹³³Xe (5.2 days) and a small leakage through mountain cracks. However, a direct number of the fission inventory is hard to estimate as the fission yield depends not only on the explosive yield but also on the type of nuclear device. Within the available data it is not possible to determine if nuclear fusion was involved in the process. The largest yield ever achieved by a boosted fission device was 720 kT TNT equivalent of the British Orange Herald in 1957. For this explosion an amount of 118 kg ²³⁵U was mounted (Arnold and Pyne, 2001). The explosive yield in this study of between 160 and 400 kT TNT equivalent is therefore still fully compatible with a fission-only device.

DInSAR results show a large area of strong displacement in the case of the 2017 test, and it is not possible to demark a definite surface location of the test epicenter from these measurements. However, retrieving a ground truth from DInSAR measurements from the January 2016 test, which showed a clear and distinct maximum of surface displacement, can help to provide an absolute epicenter location of the 2017 test. By combining remote sensing and relative locations methods only the epicenter of the event can be estimated. A third method (depth phase modeling) is required for the complete estimation of the 2017 test hypocenter. This method is further able to provide an explanation for the large difference in m_b and M_w, which has not been observed in this order for the previous tests. Furthermore while for previous weaker tests radionuclide measurements were able to prove the nuclear character of the explosions (Ross et al., 2017), radionuclide measurements alone in September 2017 could have been misleading with regard to the source area. This examples clearly demonstrate that an integrated multi-technology and multi-methodology approach is essential for a reliable and comprehensive identification and characterization of a potential nuclear test.

7 Conclusions

Overall, the following main conclusions can be deduced from this study.

• The explosive character of the 3 September 2017 North Korean event is confirmed by MTI analysis.

• The yield of the event is estimated to be around 400 kT TNT equivalent with indications that this value should be considered as an upper limit and might be as low as 160 kT TNT equivalent. This overestimation can be explained by the enhancement of peak amplitudes of teleseismic P waves due to topography and depth phase effects.

• The estimated yield of the nuclear device is certainly smaller than the largest documented yield ever achieved by a boosted fission device and is therefore still compatible with a fission-only device.

• Strong surface deformations (±10 cm) are observed in an area of 3×4 km². Comparisons of affected surface area to the January 2016 test underline the enormous power of the 2017 nuclear test. Furthermore, multiple landslides as well as a number of aftershocks were observed in the aftermath of the test.

• The aftershock directly following the nuclear explosion has a similar depth as the test and can be characterized as a collapse of the cavity created by the test.

• Infrasound signals measured at the station I45RU can clearly be related to the 2017 event. It is the first time that a thermospheric propagation path of infrasound could be observed for an underground nuclear explosion.

• No immediate measurements of radionuclides related to the test in September were observable, but later occurrences of radionuclides are consistent with a delayed leakage from the test site in October.

• The test site might have been strongly stressed and shattered and thus rendered useless for further test activities.

The combination of the results from the different technologies and methods yields a reliable estimation of hypocenter, magnitude, explosive energy, source mechanism, and an indication of delayed leakage of ¹³³Xe from the test site. The cross-links between the different results complement each other in a consistent manner. The multi-technology and multi-methodology analysis presented in this study clearly indicates that the September 2017 North Korean event was in fact a nuclear test. Even in the phase of before its entry into force, the CTBTO verification regime has again demonstrated its readiness with respect to the recent nuclear test.