We compute dynamic stress changes induced in the magma chamber by a nuclear explosion using PyLith, “a finite-element code for dynamic and quasistatic simulations of crustal deformation, primarily earthquakes and volcanoes”\cite{Aagaard_2013a,Aagaard_2013b}.
A magma chamber has been detected beneath Baekdu from low shear wave velocity in the depth range of 10-16 km \cite{ZHANG_2002}. Based on this observation, we assume that the magma chamber is a sphere of a 3 km radius and is contained in a rectangular box representing the crust (Fig. XX) We approximate the waves propagating through the magma chamber as a pulse of plane wave. Our goal to estimate only the peak stress change justifies this simplification. The plane wave is generated by the time-dependent displacement uniformly applied on one boundary, which is given as \[u(t) = u_{0}\frac{t}{\tau}\exp\left(-\frac{t}{\tau}\right),\] where \(u\) is the boundary-normal component of displacement, \(t\) is time and \(u_{0}\) and \(\tau\) are set to be 0.24 m and 1.5 sec. These values are chosen such that the peak ground velocity (PGV) becomes 0.004 m/s at the moment when the plane wave first reaches the magma chamber (Fig. XX). This value of PGV is acquired from the PGV attenuation equation constructed from the attenuation pattern around Baekdu. The opposite side of the box has an absorbing boundary condition while the sides parallel to the wave’s propagation direction are assumed to be reflecting boundaries (i.e., the boundary-normal displacement component is zero.) The stresses are initially zero everywhere in the domain and the gravitational body force is not included. These conditions are suitable to this study’s purpose of estimating only dynamic stress changes.
To represent various possible states of the magma chamber during an explosive volcanic eruption \cite[e.g.,][]{Hong_2014}, we assigned Vp/Vs ratios between 1.65 and 1.95. Vp, Vs, Vp/Vs as well as bulk and shear modulus (K and G) and Poisson’s ratio (\(\nu\)) are listed in Table \ref{tab:plane_wave_magma_chamber_parameters}. The uniform density of the magma chamber is fixed at 2500 kg/m\(^{3}\).
\label{tab:plane_wave_magma_chamber_parameters}
| Model | Vp/Vs | Vp (km/s) | Vs (km/s) | K (GPa) | G (GPa) | \(\nu\) |
| 1 | 1.65 | 5.90 | 3.58 | 44.3 | 32.0 | 0.209 |
| 2 | 1.75 | 5.90 | 3.37 | 49.2 | 28.4 | 0.258 |
| 3 | 1.85 | 5.90 | 3.19 | 53.1 | 25.4 | 0.293 |
| 4 | 1.95 | 5.90 | 3.03 | 56.4 | 23.0 | 0.321 |
The surrounding crust has Vp and Vs equal to 6.17 km/s and 3.57 km/s. The crustal density is 2700 kg/m\(^{3}\).
The process of wave generation and propagation in the model with Vp/Vs equal to 1.85 is shown in Fig. \ref{fig:plane_wave_walkthrough}.
Dynamic stress changes in the magma chamber with Vp/Vs equal to 1.85 is traced in terms of the maximums of \(\sigma_{xx}\) and pressure within the chamber, which are denoted as \(\sigma_{xx,max}\) and \(p_{max}\) (Fig. \ref{fig:plane_wave_stress_change}). The peak values of the mean and the maximum dynamic pressure are 1.48 and 1.70 kPa and are reached at about 32 seconds since the initiation of the simulation, which corresponds to the arrival of surface waves(?). Post-peak fluctuation of the dynamic pressure change is persistent as expected from the lack of attenuation mechanism in the model.
Peak values of the dynamic pressure change are positively correlated with Vp/Vs ratios (Fig. \ref{fig:pressure_vs_VpVs}). The maximum increases by about 0.25 kPa from 1.65 kPa at Vp/Vs = 1.6 to 1.90 KPa at Vp/Vs = 2.0. Increase in the peak mean value is about 0.15 KPa: 1.45 kPa at Vp/Vs = 1.6 to 1.6 kPa at Vp/Vs = 2.0.
We compute pressure changes induced in the magma chamber by a nuclear explosion using PyLith, “a finite-element code for dynamic and quasistatic simulations of crustal deformation, primarily earthquakes and volcanoes”\cite{Aagaard_2013a,Aagaard_2013b}.
The region around Baekdu and the nuclear test site is represented by a 280 \(\times\) 300 \(\times\) 100 km domain (Fig. \ref{fig:model_setup}a). The domain has a two-layer structure composed of a crustal layer with density of 2700 kg/m\(^{3}\) and a mantle layer of 3300 kg/m\(^{3}\). The crust is about 33 km thick on average in the region but gets as thick as 40 km around Baekdu \cite{Hong2008}. We adopt seismic wave speeds inferred for the region by \citet{DUAN_2005} and \citet{ZHANG_2002}: Vs = 3.25 km/s and Vp = 5.62 km/s in crust; 4.0 and 8.0 km/s in mantle.
A magma chamber has been detected beneath Baekdu from low shear wave velocity in the depth range of 10-16 km \cite{ZHANG_2002}. Based on this observation, we assume that the magma chamber is a sphere of a 3 km radius, of which center is at 13 km depth and 120 km away to the northeast of the explosion site at the center of the domain. To represent various possible states of the magma chamber during an explosive volcanic eruption \cite[e.g.,][]{Hong_2014}, we assigned Vp/Vs ratios between 1.6 and 2.0. Vp, Vs, Vp/Vs ratio as well as bulk and shear modulus (K and G) and Poisson’s ratio (\(\nu\)) are listed in Table \ref{tab:magma_chamber_parameters}. The density of the magma chamber is assumed to be uniformly 2500 kg/m\(^{3}\).
\label{tab:magma_chamber_parameters}
| Model | Vp/Vs | Vp (km/s) | Vs (km/s) | K (GPa) | G (GPa) | \(\nu\) |
| 1 | 1.600 | 5.200 | 3.300 | 31.30 | 27.23 | 0.1628 |
| 2 | 1.700 | 5.100 | 3.000 | 35.03 | 22.50 | 0.2354 |
| 3 | 1.800 | 5.040 | 2.800 | 37.37 | 19.60 | 0.2768 |
| 4 | 1.900 | 5.035 | 2.650 | 39.97 | 17.56 | 0.3084 |
| 5 | 2.000 | 5.000 | 2.500 | 41.67 | 15.63 | 0.3333 |
The above described geometries are discretized into an unstructured mesh of tetrahedra (Fig. \ref{fig:model_setup}b). Element edge lengths are about 10 km but as small as 0.2 km around the explosion site and the magma chamber (Fig. \ref{fig:model_setup}c).
The stresses are initially zero everywhere in the domain and the gravitational body force is not considered. These conditions are suitable to this study’s purpose of estimating only dynamic stress changes. Absorbing boundary conditions are applied to all the boundaries except the top boundary, which is a free surface, and the boundary for explosion.
A nuclear explosion is approximated by time-dependent normal tractions applied on the surface of a 1 km-radius sphere carved out of the crust(Fig. \ref{fig:model_setup}c). The hollow sphere is at the center of the domain and buried at the depth of 2 km. The actual explosion must have occurred at a shallower depth but we used this geometry to prevent the mesh around the explosion site from being refined too much. The magnitude of normal traction is spatially uniform over the spherical boundary but varies in time. The magnitude is 0 Pa from 0 to 0.1 sec, linearly increases to 75 MPa from 0.1 to 2.1 sec and linearly decreases to 0 Pa from 2.1 to 4.1 sec. It remains at 0 Pa afterwards. The maximum magnitude and time variation of the tractions are adjusted such that model outputs can reasonably fit waveforms recorded at the two nearby seismic stations, DNH and YNB, in terms of the magnitude and duration. The monitoring points are set at the locations corresponding to the relative coordinates of these two stations with respect to the explosion site (Fig. \ref{fig:model_setup}b). DNH is at about 210 km NNW and YNB at about 160 km NNE from the test site. Displacement solutions interpolated at the monitoring points are compared with the records of the nuclear test-generated ground motions (Fig. \ref{fig:displ_mag}).
The process of wave generation and propagation in the model with Vp/Vs equal to 1.7 is shown in Fig. \ref{fig:walkthrough}. Iso-magnitude surfaces for the displacement field reflects the main features included in the models such as the isotropic source, the layered structure composed of crust and mantle and the absorbing boundary conditions. As in the displacement magnitude plot (Fig. \ref{fig:displ_mag}), the dynamic pressure change is seen to attenuate slowly
Dynamic pressure change in the magma chamber, of which Vp/Vs is 1.7, is traced in terms of the mean and the maximum pressure in the chamber (Fig. \ref{fig:P_change_mean_max}). The peak values of the mean and the maximum dynamic pressure are 1.48 and 1.70 kPa and are reached at about 32 seconds since the initiation of the simulation, which corresponds to the arrival of surface waves(?). Post-peak fluctuation of the dynamic pressure change is persistent as expected from the lack of attenuation mechanism in the model.
Peak values of the dynamic pressure change are positively correlated with Vp/Vs ratios (Fig. \ref{fig:pressure_vs_VpVs}). The maximum increases by about 0.25 kPa from 1.65 kPa at Vp/Vs = 1.6 to 1.90 KPa at Vp/Vs = 2.0. Increase in the peak mean value is about 0.15 KPa: 1.45 kPa at Vp/Vs = 1.6 to 1.6 kPa at Vp/Vs = 2.0.