## Abstract

The 3 September 2017 $Mw$ 5.2 North Korean underground nuclear test (DPRK2017) is the largest man‐made explosion with surface displacements observed by Synthetic Aperture Radar (SAR) and showed as much as 3.5 m of horizontal permanent deformation. Although regional distance waveform‐based seismic moment tensor (MT) inversion methods successfully identify this event as an explosion, the inverted solutions do not fit the SAR displacement field well. To better constrain the source, we developed an MT source‐type inversion method that incorporates surface ground deformation (accounting for free‐surface topography), regional seismic waveforms, and first‐motion polarities. We applied the source‐type inversion over a grid of possible source locations to find the best‐fitting location, depth, and point‐source MT for the event. Our best‐fitting MT solution achieves ∼70% horizontal geodetic fit, ∼80% waveform fit, and 100% fit in the first‐motion polarities. The joint inversion narrows the range of acceptable source types improving discrimination, and reduces the uncertainty in scalar moment and estimated yield. The method is transportable and can be applied to other types of events that may have measurable geodetic signals such as underground mine collapses and volcanic events.

## Introduction

On 3 September 2017, the Democratic People’s Republic of Korea conducted its largest nuclear test (DPRK2017) at the Pungee‐Ryi test site having an $mb$ 6.3 (U.S. Geological Survey [USGS] and National Earthquake Information Center) and a scalar moment of $8.06\xd71016\u2009\u2009N\xb7m$ corresponding to $Mw$ 5.2 (Chiang *et al.*, 2018). The event was located beneath Mt. Mantap, which geologically is inferred to be granodiorite at emplacement depth that is nonconformably overlain by stratified volcanic rocks with a basaltic cap (Pabian and Coblentz, 2018). Figure 1a shows the region and locations of seismic stations.

This event had a substantially larger yield than previous tests and provides a high‐quality seismic waveform dataset from the Chinese Digital Seismic Network (BJT, HIA, and MDJ), Korean Seismic Network (TJN), and Incorporated Research Institutions for Seismology (IRIS)‐USGS Global Seismographic Network (GSN; INCN and MAJO) stations. These data have been utilized in studies of DPRK2017 and previous explosions at the test site (e.g., Ford *et al.*, 2009, 2012; Alvizuri and Tape, 2018; Chiang *et al.*, 2018). Above‐shot ground displacement obtained from pixel tracking of Synthetic Aperture Radar (SAR) data from the German TerraSAR‐X satellite (Wang *et al.*, 2018; Fig. 1b) reveals as much as 3.5 m of horizontal displacement (generally radially outward from the epicenter), and a heterogeneous pattern of localized relatively low‐amplitude uplift and subsidence (∼0.5 m) at the test site. Wang *et al.* (2018) inverted the ground deformation observations by finding a model consisting of an explosion at 450 m depth, an associated collapse, and material bulking in the rock volume surrounding the shot point. They located the epicenter at the crest of Mt. Mantap, consistent with travel‐time locations (e.g., Gibbons *et al.*, 2018; Myers *et al.*, 2018). There was a seismic collapse event 8.5 min after the explosion, which had a moment less than 1/10 of the explosion moment (e.g., Chiang *et al.*, 2018), and possible surface deformation associated with it is negligible compared to that from the main explosion. Hydrodynamic simulations of Stevens and O’Brien (2018) illustrate reduced vertical motions due to topographic effect, gravity, and nonlinear processes such as spall and cracking.

Chiang *et al.* (2018) modeled the regional waveforms and first‐motion polarities by finding the seismic moment tensor (MT) and identifying the event as an explosion. However, as shown in Figure 1c, their solution does not fit the SAR data well (The details of the calculations are described in the following sections.). Dreger *et al.* (2018) found a subset of the acceptable long‐period waveform, first‐motion source‐type inversion solutions that can better explain the above‐shot displacements, motivating the development of the joint inversion of SAR static ground deformation, regional seismic waveforms, and first‐motion polarities presented here.

Such joint inversions have become commonplace in the study of large earthquakes (e.g., Wald *et al.*, 1996), and it has been shown that the seismic waveform and geodetic datasets are complementary in providing overall greater constraint on the source process. We develop a joint inversion scheme for the DPRK2017 dataset, modifying the source‐type‐specific MT inversion method (Nayak and Dreger, 2015) incorporating the surface displacements in addition to the regional distance waveforms and first‐motion polarities. We use static displacement Green’s functions that account for significant topography at the test site. Although there is source nonlinearity (e.g., Stevens and O’Brien, 2018), the horizontal surface deformation is generally consistent with the seismic observations and the model of an elastic seismic source (Fig. 1; Dreger *et al.*, 2018). Therefore, it is worth evaluating whether an elastic‐model joint inversion, which can be applied rapidly after the event, might improve recovery of source parameters and improve source‐type discrimination capability. As follows, the joint inversion does find an MT solution that satisfies the three datasets and results in a narrower range of acceptable solutions improving discrimination as well as reducing uncertainty. The developed method is transportable to other areas of interest in which substantial near‐source ground deformation may be observed by SAR satellites and may also be applied to other types of problems that involve complex source mechanics such as volcanic eruptions.

## Inversion Method

To solve for the MT of the event, we modify the source‐type‐specific inversion method of Nayak and Dreger (2015), which finds the best‐fitting eigenvectors and a moment scale factor for specific ratios of eigenvalues (source‐type) uniformly distributed over source‐type space (e.g., Hudson *et al.*, 1989; Tape and Tape, 2012). The method essentially maps the goodness of fit to the source‐type solution space, enabling improved discrimination and assessment of solution uncertainties (e.g., Ford *et al.*, 2010; Nayak and Dreger, 2015; Alvizuri and Tape, 2018). Because of the dominance of surface waves in regional distance waveforms for shallow focus events and an azimuthally uniform source radiation pattern for explosions, there is a trade‐off between pure explosion source types and vertically oriented compensated linear vector dipoles (CLVDs) with a major vertical compression axis (Ford *et al.*, 2010). This is remedied by jointly inverting regional waveforms and regional and teleseismic first‐motion polarities (Ford *et al.*, 2012; Chiang *et al.*, 2014). In this study, we incorporate the static deformation field from SAR observations as a third dataset in a modified joint inversion based on the method of Nayak and Dreger (2015) to further constrain the source mechanism of DPRK 2017. We used the SAR pixel‐mapping ground deformation observations from Wang *et al.* (2018). Inversions are performed at 672 possible source locations (Fig. 2a) sampled 200 m horizontally and 100 m vertically to find the best‐fit location for the deformation data. This grid is centered on the epicenter from Myers *et al.* (2018) with 2 km × 2 km horizontal dimension and 900 m vertical extent that encompasses the locations of the second through sixth DPRK nuclear tests and covers the full range of possible source depth limited by tunnel elevation.

The waveform data were instrument corrected, integrated to ground displacement, and band‐pass filtered (50–20 s, except MAJO 50–33 s period) using an acausal four‐pole Butterworth filter. Green’s functions for the regional distance waveforms were computed using the MDJ2 velocity model (Ford *et al.*, 2009) using the Herrmann (2013) frequency–wavenumber integration (FKI) code. We picked a total of 50 *P*‐wave polarities on vertical‐component regional and teleseismic waveforms from IRIS (GSN and China Digital Seismograph Network) and F‐Net stations. The takeoff angles for modeling the polarities were calculated using the IASP91 velocity model.

Surface waves can be affected by topography, and, as shown by Stevens and O’Brien (2018), this can partially explain the $Ms$ anomaly observed for DPRK nuclear tests in which $Ms$ is larger than for flat‐surface cases. Nevertheless, it has been shown that surface waves are reasonably modeled in the 0.02–0.05 Hz passband in MT inversions using plane‐layered velocity structure and are able to discriminate events as explosions (Ford *et al.*, 2009; Alvizuri and Tape, 2018; Chiang *et al.*, 2018). Although beyond the scope of this study, future work may consider utilizing a common 3D model for the geodetic and regional waveform Green’s functions to better address this issue.

We compute the static deformation Green’s functions using the SW4 summation‐by‐parts finite‐difference seismic‐wave propagation code that allows for free‐surface topography (Rodgers *et al.*, 2010; Petersson and Sjögreen, 2012, 2015; Sjögreen and Petersson, 2012). The elastic model for the SW4 calculations is defined as a homogeneous half‐space ($VP=2.6\u2009\u2009km/s$, $VS=1.5\u2009\u2009km/s$, and $\rho =2570\u2009\u2009kg/m3$). Topography was incorporated from the Shuttle Radar Topography Model with a 30 m sampling. This model has slower velocities than MDJ2, which gives the regional average velocity values and is based on initial modeling of the magnitude of the static deformation field using regional waveform MT models as a constraint, accounting for the combined effect of the local site conditions of weathered granodiorite and stratified volcanic rocks at the test site. Green’s functions were computed at the 672 possible source points individually for each of the six MT elements. The finite‐difference grid spacing was 25 m with a 10 × 10 × 6 km computational domain centered on the event location (Myers *et al.*, 2018), comprising 49.2 million grid points. We computed 16 s seismogram durations for each location and MT component. Each calculation ran in about 1 min using four nodes of the Lawrence Livermore National Laboratory Lassen graphic processing units platform. The static displacements were determined as the average of the final 5% of the simulated time series. We verified flat half‐space model static displacements calculated with SW4 by comparing with those calculated using FKI (Herrmann, 2013).

## Results

The joint inversion results for the 672 possible source locations achieves a maximum of about a 70% variance reduction in horizontal‐component geodetic fit shown in Figure 2b,c. The results in Figure 2 were determined using the joint inversion with fixed weights of 40% waveform weight, 10% first‐motion weight, and 50% geodetic weight. The epicentral location we obtain is the same as reported in other studies within uncertainty (Myers *et al.*, 2018; Wang *et al.*, 2018); however, the depth of the event is approximately 400 m shallower than Myers *et al.* (2018) and outside their uncertainty estimate. The solution shown Figure 2 is close to that found by Wang *et al.* (2018). The shallower depth, which better satisfies the geodetic data, may be due to finite‐source effects of the explosion and vertically asymmetrical damage (e.g., Stevens and O’Brien, 2018) to which the surrounding medium responds, as well as other processes such as chimneying, spallation, and related deformation (e.g., Patton and Taylor, 2011). We discuss the effect of source depth on the inversion results in the Discussion section.

The relative weights among the three datasets in the joint inversion are free parameters that need to be evaluated. To do this, we performed 14,553 inversions (63 source locations; see Fig. 2a, each with 231 relative data weight combinations.). These selected sources are the closest 63 locations from the epicenter reported by Myers *et al*. (2018). The weight of each data type varied between 0% and 100% of the summed weight of the three datasets in intervals of 5%. Figure 3a shows the results of those inversions plotted against the horizontal geodetic and three‐component waveform fit. Individual trade‐off trends can be seen for the different possible source positions. The red‐dashed lines show 95% and 90% fit thresholds for the geodetic and waveform data, respectively, identifying the subset of best‐fitting joint solutions. Figure 3b provides a zoomed‐in view of this subset of 33 solutions, 32 of which are for the same source point shown in Figure 2. All of these solutions fit the first motions perfectly with 100% variance reduction, which is consistent with the fact that first motions are mostly unaffected by nonlinearities and topography.

There are two main clusters of solutions identified. The four solutions located in the upper part of Figure 3b correspond to MTs with larger isotropic components that improve the horizontal geodetic fit by 1%–2% at the expense of vertical geodetic fit. It is likely that these two populations represent local minima in the linearized source‐type inversion. To examine this, we ran the inversion 100 times each time with three starting models giving a total number of 300 different randomized initial solutions to test for convergence. As Figure 3c shows, many of the solutions in the lower population rise to the upper cluster. For both clusters, the range in waveform and first motion fit is similar. We chose a solution from the upper cluster because of the better horizontal geodetic fit.

Figure 4 shows the distributions of the six MT components for all of the solutions in the upper cluster of Figure 3c. Except for the $Mxy$ component, all components have an uncertainty of less than 4%, represented by two standard deviations of the distribution. Statistically, the solutions of the upper cluster show little variation, and we chose our preferred solution based on the minimum Euclidean distance from the average of all the inversions in the upper cluster. Figure 5 shows the fit of the preferred solution, which has 79.2% waveform and 69.9% horizontal geodetic variance reductions. This solution has a dominant isotropic component (66.3% of the total seismic moment), and 9.6% double couple (DC), and 24.7% CLVD using the decomposition, assuming the same principal stress orientations for the DC and CLVD components (Jost and Herrmann, 1989).

## Discussion

The preferred solution has a very good waveform fit (Fig. 5a), especially for the radial and vertical components as expected for an explosion. The slightly inclined compressional CLVD fits the Love waves well, including the tangential component radiation nodes at stations HIA and MAJO.

The fit to the horizontal deformation (Fig. 5b) is also very good. As expected for an explosion, the deformation is radially outward at all azimuths, with the largest offsets toward the west and south on the steep slopes of Mt. Mantap. There is a strong azimuthal variation in the amplitude of the horizontal deformation that is due to steep western and southern slopes of Mt. Mantap. Note the substantially smaller observed and simulated motions at similar distance sites to the northeast. We compared flat surface and topographic Green’s functions, and found that displacements correlated with topography in which larger horizontal displacements occur at sites on the steep slopes. Though the horizontal geodetic fit is about 70%, the total geodetic fit (including the vertical displacement) is 48%. We are unable to explain the heterogeneous low‐amplitude vertical deformation entirely and are unsuccessful in fitting the subsidence directly above the shot point, which, as shown by Stevens and O’Brien (2018), could be due to nonlinearity and the effects of gravity, as well as bulk compaction after the explosion (e.g., Wang *et al.*, 2018), which may have taken place in the days following the explosion (The latest scene was for 12 September, nine days after the explosion.). The erasure of the expected ∼3.5 m of uplift in the observed deformation suggests (e.g., fig. 1c of the Chiang *et al.*, 2018 MT solution) that the vertical component may be more strongly affected by complexities such as spall, chimneying, cavity collapse (Although estimated deformation for the collapse event is less than 1/10 for the explosion.), and other nonlinear processes. We find that the $Mzz$ component is smaller than either $Mxx$ or $Myy$, with the decomposition indicating a dominant explosion and approximately ⅓ of the moment from a compressional CLVD, which could represent spall‐related compaction. Although the actual source involves nonlinear processes, such as bulk dilation, vaporization, melting, and deformation of the cavity due to ambient stress, the mechanism found from only the near‐field deformation data is consistent with the far‐field observations and when jointly inverted provides additional constraint on the overall source mechanism.

The recovered depth is significantly shallower than found from travel times (Myers *et al.*, 2018). This is primarily an effect of the geodetic data, because the regional long‐period waveforms are relatively insensitive to source depth. It is possible that considering layered velocity structure within the test site might change the depth estimates and will be the focus of future work. The most likely emplacement depth is consistent with the adit elevation plus a small (1%–2%) tunnel slope for drainage giving about 730 m of overburden. Myers *et al.* (2018) find a depth of 760 ± 200 m. We examined models within their uncertainty range, and found similar source‐type plots and solutions. For example, a source with a depth of 550 m results in a solution only 7% larger in total scalar moment and 8% larger in isotropic moment, both within the estimated uncertainties of our preferred solution but with a 5% reduction in fit to the geodetic data. The $Mzz$ component does tend to increase with increasing source depth due to better fitting the three‐component SAR data. Interestingly, this trend is opposite to the effect of vanishing traction on MT recovery found by Chiang *et al.* (2016), suggesting the joint inversion can help to reduce the effects of vanishing traction on the long‐period regional Green’s functions.

The joint inversion greatly improves the discrimination of the event, as shown in Figure 5. There is a much more tightly focused region of acceptable source‐type solutions closer to the pure‐explosion pole, which is also true for the other source depths at the epicentral location. The joint inversion results in an overall lower scalar moment, and estimated yield are compared to the waveform and first‐motion joint inversion. Considering the solutions that are within 95% of the respective maximum fit, we find that the uncertainty in the isotropic and Bowers and Hudson (1999) total scalar moments as well as the estimated yield is significantly reduced (Fig. 5c,d). Using the isotropic moment, its uncertainty, an effective source depth of 500m, the MDJ2 elastic model, and an assumed gas porosity of 1%, we estimate the yield to be 191.8 ± 46.9 kT using the Denny and Johnson (1991) model. Walter and Wen (2018) reported that the range of published yield is from about 100 to 300 kT. Our estimate including 2‐sigma uncertainty is less than this overall range.

## Conclusions

A joint source‐type inversion utilizing SAR deformation data, regional seismic waveforms, and first‐motion polarity was applied to the DPRK2017 nuclear test. It was found that the joint inversion resulted in a better resolved dominantly explosive source‐type, and reduced uncertainty in scalar moment and estimated yield. This analysis shows that the static displacement field can reduce the range of possible solutions as represented in source‐type space and therefore improve MT‐based discrimination capability for large yield explosions. Although the method assumes linear elasticity, it can be performed relatively rapidly, delayed only by postevent satellite data recovery enabling timely source‐type discrimination capability. Future work at the Punggye‐ri test site should consider layered elastic models incorporating what is known about geological composition of the test site (e.g., Coblentz and Pabian, 2015; Pabian and Coblentz, 2018) to see if improvements in source depth can be obtained. Hydrodynamic modeling with nonlinear geomechanics using the source reported here could be investigated to better understand near‐source processes and deformation (Xu *et al.*, 2014; Stevens and O’Brien, 2018).

This joint inversion method shows great promise in better constraining the MT of shallow buried explosions, and because it uses satellite‐based topography and ground deformation, it is transportable to other regions with calibrated regional and local elastic models. The combination of regional seismic waveforms and near‐field deformation data is complementary and provides additional constraint on the source mechanism of the event. Although we applied this approach to the DPRK2017 nuclear test, it may also be used in other applications such as underground mine collapses and volcanic events in which there may be measurable static displacements.

## Data and Resources

Seismic waveform data used in this study were downloaded from the Incorporated Research Institutions for Seismology (www.iris.edu). The Synthetic Aperture Radar (SAR) deformation data from 3D pixel mapping is available online as a supplement to Wang *et al.* (2018). SW4 is an open‐source anelastic finite‐difference code for simulating seismic motions and is available at https://geodynamics.org/cig/software/sw4/. All websites were last accessed in August 2021.

## Declaration of Competing Interests

The authors acknowledge that there are no conflicts of interest recorded.

## Acknowledgments

The authors thank Carl Tape, Jeff Stevens, and Steven Gibbons for their constructive comments and suggestions that improved this article. Rodrigo Chi‐Durán acknowledges the support from the National Agency for Research and Development (ANID) Scholarship Program/DOCTORADO BECAS CHILE/2015-2305615000. RCD also thanks Leonel Sanchez for his advice about figures. The authors are grateful to the Livermore Computing Grand Challenge Program for access to the Lassen graphic processing units (GPU) platform for SW4 simulations. This work was supported by the Department of State under Contract Number 19AQMM19P1591. This work was performed in part (Arthur J. Rodgers) under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract Number DE-AC52-07NA27344.