Prediction of ground motion and dynamic stress change in Baekdusan (Changbaishan) volcano caused by a North Korean nuclear explosion

Strong ground motions induce large dynamic stress changes that may disturb the magma chamber of a volcano, thus accelerating the volcanic activity. An underground nuclear explosion test near an active volcano constitutes a direct treat to the volcano. This study examined the dynamic stress changes of the magma chamber of Baekdusan (Changbaishan) that can be induced by hypothetical North Korean nuclear explosions. Seismic waveforms for hypothetical underground nuclear explosions at North Korean test site were calculated by using an empirical Green’s function approach based on a source-spectral model of a nuclear explosion; such a technique is efficient for regions containing poorly constrained velocity structures. The peak ground motions around the volcano were estimated from empirical strong-motion attenuation curves. A hypothetical M7.0 North Korean underground nuclear explosion may produce peak ground accelerations of 0.1684 m/s2 in the horizontal direction and 0.0917 m/s2 in the vertical direction around the volcano, inducing peak dynamic stress change of 67 kPa on the volcano surface and ~120 kPa in the spherical magma chamber. North Korean underground nuclear explosions with magnitudes of 5.0–7.6 may induce overpressure in the magma chamber of several tens to hundreds of kilopascals.

Various studies have reported on volcanic eruptions triggered by earthquakes [21][22][23][24][25] . There is growing anxiety whether a large future nuclear explosion at the North Korean test site may disturb the magma chamber and cause a volcanic eruption. Baekdusan is located 116 km away from the North Korean test site. It is close enough to be potentially affected by a moderate-sized or large seismic event 26 . The triggering of a volcanic eruption is controlled by various parameters including a stress (pressure) change, the magma composition, the depth and volume of the magma chamber, the viscosity, and local tectonic settings 25,27 . The stress (pressure) change plays a crucial role in volcanic eruption.
The internal pressure in a magma chamber can be influenced by a dynamic stress change as well as lithostatic stress field 28 . Previous UNE tests have shown that detonation strength can be greater than M7.0 29 . Large UNEs may produce strong seismic waves over local and regional distances. Dynamic stress changes induced by transient seismic waves can overpressurize the magma and trigger volcanic eruption [30][31][32][33] .
We estimate dynamic stress changes that a hypothetical North Korean UNE can induce in the Baekdusan magma chamber. However, the Korean Peninsula on which Baekdusan is located was formed by a series of continental collision and rifting from the late Permian to mid-Miocene [34][35][36] , which led to the construction of complex crustal structures 17,34,37 (see, supplementary materials). Synthetic calculation of seismic waveforms is not suitable for such regions where the velocity structures are poorly constrained since seismic waveforms are sensitive to the crustal structure and surface topography along the raypath 13,19,38 . We calculate the seismic ground motions from hypothetical nuclear explosions by using an empirical Green's function approach based on the seismic waveforms of previous nuclear explosions and a UNE source-spectral model. The peak ground velocities and dynamic pressure changes in the magma chamber of Baekdusan volcano are determined from a numerical model of the dynamic stress coupling and distance-dependent ground motion attenuation.

UNE source calibration.
A UNE source-spectral model [39][40][41] is implemented to calculate the seismic waveforms for a given UNE detonation size. The validity of the UNE source-spectral model is checked, and the apparent seismic moment in the model is calibrated for the UNE magnitude 20 . The spectral-amplitude ratios of seismic records of North Korean UNEs at common stations yield source-spectral amplitude ratios of the UNEs because raypath effects such as geometrical spreading and seismic attenuation are naturally corrected.
The spectral-amplitude ratios between the 2009 and 2013 UNEs are calculated for six azimuthal ranges, 70-95°, 95-120°, 120-145°, and 145-178° for stations on Japanese islands, and 176-200° and 200-227° for stations on the Korean Peninsula ( Fig. 1(a-c)). The spectral-amplitude ratios are found to be similar for all azimuths ( Fig. 1(c)). The comparison between the observed and theoretical spectral-amplitude ratios verifies the theoretical UNE source-spectral model. The overshoot parameter in the source-spectral model is ξ = 1.05 20 . The seismic waveforms are corrected for instrumental responses and bandpass-filtered between 0.01 and 30 Hz. The zero-to-peak ground motions (accelerations, velocities) are measured. A reference strong-motion attenuation curve is determined as a function of distance from the observed peak ground motions. The event-strength-calibration constants adjusting the levels of peak ground motions for hypothetical UNEs are determined from the quasi-observed waveforms ( Peak dynamic stress change in magma chamber. We compute the dynamic stress changes induced in the magma chamber by an m b 7.0 nuclear explosion with PyLith, which is a finite-element code for dynamic and quasi-static simulations of crustal deformation 42 . We consider an impulsive plane wave incident to the magma chamber in order to estimate the upper bound of the peak dynamic stress change. The plane wave is generated in the left margin of the domain. The strength of the input signal is set to produce transient waves of which peak ground motion is equivalent to that expected from quasi-observed seismic waveforms on the volcano surface (Fig. 4).
A spherical model with a radius of 3 km and three spheroidal models with a flattening of 0.5 are considered for the magma chamber ( Fig. 4(a)). We assess the variation in the induced stress depending on the geometry and density of the magma chamber. The densities of the magma chamber are set to be 2500, 2580 and 2660 kg/m 3 for three different compositions of the magma chamber, representing 50, 30 and 10 vol% of silicic magma compared to the typical continental crust (2700 kg/m 3 ) 43 .
To represent various possible states of the magma chamber during an explosive volcanic eruption 44 , we consider four V P /V S ratios of 1.65, 1.75, 1.85, and 1.95 (Table 1), representing various magma compositions and pore-fluid pressures. The V P /V S ratio is low in felsic rock and high in mafic rock 45 . Also, the V P /V S ratio increases with pore-fluid pressure. These V P /V S ratios are equivalent to Poisson's ratios of 0.209, 0.258, 0.293, and 0.321. The transient wavefield comprised compressional pulses satisfying the observed peak ground motions on the volcano surface to induce dynamic stress changes in the magma chamber ( Fig. 4(b-d)). The radial peak dynamic stress change in the magma chamber (Δ σ xx,max ) are found to be constant for V P /V S ratios, while the peak dynamic pressure change in the magma chamber (p max ) increase with the V P /V S ratio. On the other hand, Δ σ xx,max and p max generally increase with the density. The spherical magma chamber model produces larger Δ σ xx,max and p max than the three spheroidal models.  The observations suggest that Δ σ xx,max and p max are only significant in the model Oy, which are 46 and 83 kPa, respectively, compared to 70 and 124 kPa for the corresponding model with a spherical chamber (V P /V S = 1.75, density = 2660 kg/m 3 ). The other two models, Ox and P produce smaller Δ σ xx,max and p max than model Oy. This suggests that the spherical model constitutes an environment that produces the upper bound of stress levels by incident seismic waves. The induced peak dynamic pressure change in the spherical magma chamber model ranges between 60 and 80 kPa for various V P /V S ratios and densities observed in the crust. For the magma chamber with a density of 2500 kg/m 3 , p max is estimated to be 64 kPa for a V P /V S ratio of 1.65 and 75 kPa for a V P /V S ratio of 1.95 (Fig. 4).

Discussion and Conclusions
The strong motions from hypothetical UNEs were simulated using UNE source-spectral models and seismic records of previous UNEs to avoid the errors associated with uncertainty in crustal models. The method was verified by comparisons between observed and synthetic seismic waveforms for previous UNEs. The peak ground velocities can be converted to peak dynamic stress changes. The peak dynamic stress change was found to linearly increase with the magnitude (Fig. 3(b)). The peak dynamic stress changes on the volcano surface were expected to be 3.8-157.2 kPa in the horizontal direction and 3.4-90.0 kPa in the vertical direction for hypothetical UNEs with magnitudes of 5.0-7.6 ( Fig. 3(b-d); see supplementary materials). In particular the horizontal peak dynamic stress changes on the volcano surface reached 16 kPa with an m b 6.0 UNE and 67 kPa with an m b 7.0 UNE.
The radial peak dynamic stress change (Δ σ xx,max ) and induced peak dynamic pressure change in the magma chamber (p max ) were found to generally increase with the density. On the other hand, Δ σ xx,max was constant for different V P /V S ratios, while p max increased with the V P /V S ratio. The spherical model constituted the conditions for the upper bounds of Δ σ xx,max and p max to be reached.
The Δ σ xx,max and p max values induced by an m b 7.0 UNE in spherical magma chamber models with densities of 2580 and 2660 kg/m 3 were about 3 and 5% lower, respectively, than those in the spherical magma chamber with a density of 2500 kg/m 3 . The radial peak dynamic stress change in the magma chamber with a V P /V S ratio of 1.85 and density of 2500 kg/m 3 reached the maximum of 118 kPa with averages of around 100-110 kPa during transmission through the chamber (Fig. 4(d)). The overall radial peak dynamic stress changes were generally around 110-130 kPa. The peak dynamic pressure changes were less than 80 kPa for most magma compositions.
Earthquakes are remotely triggered by peak dynamic stress changes of 0.01-1 MPa 46 . The triggering of volcanic eruption is controlled by bubbles that nucleate at different levels of stress (pressure) changes depending on the medium texture 47 . The static pressure in the magma is on the order of several hundred megapascals [47][48][49] . Reports have suggested that bubbles typically nucleate at magma chamber overpressures of several dozens of MPa 32,43,47 . Bubble nucleation may start at a pressure change of lower than 1 MPa in a medium with a microlite texture 47 . North Korean nuclear explosions are expected to produce pressure changes of tens to hundreds of kilopascals, causing concern over the possible triggering of volcanic eruption in textured media.

Methods
Quasi-observed waveform synthesis. The inter-event distances with respect to the dominant frequencies of regional seismic waves allow us to assume UNEs to be doublet seismic events. The raypath effects to common stations can be assumed to be the same between doublet events. The UNE source spectrum S( f ) can be represented as a function of the seismic moment [39][40][41] (see, supplementary materials). The source-spectral ratio between two events at a common station is given by  These relationships are particularly effective for UNE sources because they are hardly affected by azimuth-dependent source effects such as directionality and radiation patterns. The spectra of seismic waveforms for an UNE of a certain size can be calculated by using the spectra of observed seismic waveforms for a previous UNE with (1).
Peak ground motion attenuation curve. The peak ground acceleration (PGA) and peak ground velocity (PGV) attenuate with distance, and generally satisfy the relationship 50,51   =  −  −  ,  ( )   , , ,  ,,  ,  , , , where G i,j,k,l (i = PGA, PGV, j = h, v) is the peak ground motion (PGA or PGV) for the horizontal or vertical component at station k for event l in a hypocentral distance of r k,l , A i,j,l is a calibration constant for the event size, B i,j is a calibration constant for the geometrical spreading effect, and C i,j is a calibration constant for anelastic absorption. The PGA is in m/s 2 , the PGV is in m/s 2 , and the distance r is in kilometers. The calibration constants are determined so that they yield the minimum squared error (see supplementary materials). The peak ground motion attenuation curve is applied to determine the strength of a strong ground motion at a given distance.
Estimation of dynamic stress change. The radial peak dynamic stress change σ r induced by transient seismic waves is given by 52,53 where μ is the shear modulus,  u r is the peak ground velocity, and β is the shear wave velocity. We set the shear modulus at 34.95 GPa and the shear velocity at 3.58 km/s to consider the crustal properties of the Korea Peninsula at a depth of 10 km 45,54 (see supplementary materials).
The dynamic stress changes in the magma chamber is calculated using by PyLith 42 . We consider a spherical model and three spheroidal models for the magma chamber ( Fig. 4(a)). The radius of the spherical model is set to 3 km to account for the dimensions of low velocity anomalies revealed in a seismic exploration study 9 . The spheroidal models represent possible geometric orientations of a magma chamber deviating from a symmetric shape (spherical model). The spheroidal models have a flattening of 0.5. The spheroidal models comprise one prolate spheroid (P) and two oblate spheroids for which the symmetry axes are parallel with the x (model Ox) or y (Oy) axis. We consider a 12 × 15 × 15 km rectangular domain to represent a homogeneous medium that sufficiently accommodates a magma chamber model in the center.
We design a situation in which an impulsive plane wave approaches the magma chamber in order to estimate the peak dynamic stress change. The plane wave is generated in the left-hand side of the domain. The right-hand side boundary of the domain is set to an absorbing boundary condition. The boundary-normal displacement component is set to be zero on the upper and lower boundaries. The initial stress field is set to be zero over the domain, and the gravitational body force is not applied. The magma chamber is treated as an elastic body to estimate the possible upper limit of peak stress changes. The seismic velocities of the medium surrounding the magma chamber are V P = 6.17 km/s and V S = 3.57 km/s, and the density was 2700 kg/m 3 .