Traveling Ionospheric Disturbances Triggered by the 2009 North Korean Underground Nuclear Explosion

Underground nuclear explosions (UNEs) can induce acoustic-gravity waves, which disturb the ionosphere and initiate traveling ionospheric disturbances (TIDs). In this paper, we employ a multi-step and multi-order numerical difference method with dual-frequency GPS data to detect iono-spheric disturbances triggered by the North Korean UNE on 25 May 2009. Several International GNSS Service (IGS) stations with different distances (400 to 1200 km) from the epi-center were chosen for the experiment. The results show that there are two types of disturbances in the ionospheric disturbance series: high-frequency TIDs with periods of approximately 1 to 2 min and low-frequency waves with period spec-trums of 2 to 5 min. The observed TIDs are situated around the epicenter of the UNE, and show similar features, indicating the origin of the observed disturbances is the UNE event. According to the amplitudes, periods and average propagation velocities, the high-frequency and low-frequency TIDs can be attributed to the acoustic waves in the lower iono-sphere and higher ionosphere, respectively.


Introduction
Traveling ionospheric disturbances (TIDs) are the ionospheric signatures of waves and are widespread phenomena with many different origins (Afraimovich et al., 2009;Kelley, 2011;Hernandez Pajares et al., 2012).Atmospheric acoustic waves, which can be acoustic-gravity waves if the frequencies are low enough, or gravity waves associated with geophysical events, such as earthquakes (Artru et al., 2004;Liu et al., 2010), tsunamis (Artru et al., 2005;Occhipinti et al., 2006;Rolland et al., 2010), surface explosions (Calais et al., 1998) and underground nuclear explosions (Rudenko and Uralov, 1995;Krasnov and Drobzheva, 2005) can be such origins.The atmospheric acoustic or gravity waves transfer the energy from the ground into the atmosphere, which subsequently interacts with ionospheric plasma and produces disturbances.
Underground nuclear explosions (UNEs), as a disruptive event on the Earth's surface, can induce TIDs, which have gained substantial attention.Krasnov and Drobzheva (2005) present a model to describe the coupling process between radiated acoustic fields from the spall zone of UNEs and the ionosphere.The model results are consistent with highfrequency Doppler sounding data measured following underground nuclear explosions at the Semipalatinsk Test Site.Owing to notable advantages, GPS measurement is a popular technique for ionospheric sounding, including disturbance detection.Park et al. (2011) and Yang et al. (2012) used GPS observations to study the TID propagation characteristics of North Korean UNEs.
On 25 May 2009, North Korea conducted its second UNE at the reported location 41 • 17 38 N, 129 • 04 54 E. TIDs were observed after the event.Park et al. (2011) applied the numerical difference method to extract the TID series using GPS observations from five IGS stations and six stations in South Korea's local GNSS network.This method is easy to process and can magnify the amplitude of TID waves, which is efficient for weaker signatures.However, the influence of the time step on that method is not concerned.The time step could be relevant to the detectable period of spectrum.Consequently, TIDs whose periods are outside the spectrum may be "neglected".Yang et al. (2012) combined wavelet analysis and band-pass filters to obtain TID series with GEONET GPS data.This method can be used to identify TIDs with different ranges of period, but it is not easy to determine the Published by Copernicus Publications on behalf of the European Geosciences Union.
X. Zhang and L. Tang: Traveling ionospheric disturbances proper filter that can effectively distinguish the TID wave of interest from group waves, such as signal trend and noise.
In this study, we present a modified numerical difference method with multiple time steps by using GPS dualfrequency data to detect TIDs induced by the 2009 North Korean UNE.

TID detection method
As is well-known, the geometry-free combination of the dual-frequency GPS carrier phases can be used to calculate the ionospheric total electron content (TEC): where L 4 is the geometry-free combination (in meters); s represents the ionospheric TEC (in the unit TECU); k is the conversion factor between TEC and observation (k ≈ 0.105 m TECU −1 ); and b is an unknown constant bias for each satellite-receiver continuous arc of data, including carrier phase ambiguity and hardware delays.Although Eq. ( 1) cannot acquire the absolute value of TEC at a particular time, it can used to calculate the TEC variation over time with high precision, which is the basis for TID detection.Hernandez-Pajares et al. (2006) and Park et al. (2011) applied a numerical difference method to extract TIDs with order of 1 and 3, respectively.Generally, the n-order numerical difference procedure is defined: where n s(t) is n-order difference of TEC series; t is the observation epoch and τ is the time step.For a given TID with period of T , a complex exponential function can be used to express the wave d(t): where A 0 is the TID amplitude and j = √ −1 is the complex unit.Assuming the TEC trend approximately varies linearly in a short time (Hernandez-Pajares et al., 2006) and is eliminated in Eq. ( 2), then the expression between d(t) and n s(t) is obtained as follows: n s(t) = 2 n sin 2n (π τ/T )d(t). (4) According to Eq. (4), we can see the observed disturbance has the same period of d(t).That means we can pick up the disturbance waves by using Eq.(4).The augmentation factor is 2 n sin 2n (π τ/T ), which is the function of T and τ .Setting n to 2, the relationship between amplification factor and T is shown in Fig. 1 for a given τ .As seen from Fig. 1, the wave amplitude is magnified (the maximum magnification is up to 4) with an amplification factor larger than 1, which is very effective for TID detection, and vice versa.Obviously, the detectable span of disturbance period is closely associated with the time step τ and is generally limited for a specific single time step.For example, a span with amplification factor larger than 1 is only 40 to 120 s (this is the same when the order is set to 3) for τ = 30 s, which may miss other important TIDs.
To solve this problem, we use multiple time steps to replace the single time step to extend the scope of the detectable period.According to Fig. 1, the ranges for amplification factors larger than 1 where τ = 30, τ = 90 and τ = 300 s are 0.7-2, 2-6 and 6.6-20.0min, respectively.Therefore, the combined range is 0.7-20.0min (the amplification factor between spans 6 to 6.6 min is smaller than 1 but larger than 0.8, which is also effective) after applying the multi-step process.In the next section, we will use the multi-step multi-order numerical difference method to detect the TIDs induced by the North Korean UNE on 25 May 2009.As the 2-order difference is sufficiently sensitive to weaker signals, and 3order difference may lead to significant ambient noise when a longer step is used, such as 300 s, the magnitude of difference order is set to 2.

Observations and results
Eight IGS GPS stations located around the epicenter are selected to observe the ionospheric disturbances induced by the North Korean UNE on 25 May 2009 (see Fig. 2).The data sampling rate of GPS observations is 30 s. Stations CHAN and DAEJ are near the epicenter, with distances of 400-500 km, and the other stations are farther away, with distances of 900-1200 km.The height of the single-layer ionosphere is set to 350 km (the altitude of maximum ionospheric electronic density on the event day, obtained by International Reference Ionospheric model) to calculate the ionospheric TEC and location of pierce points.According to the Space Weather Prediction Center (SWPC, http://www.swpc.noaa.gov/), the Kp index was relatively low (< 1), which implies low geomagnetic activity during the time of the UNE.
Figure 3 shows the ionospheric slant TEC (STEC, which contains constant bias) calculated using the dual-frequency GPS data of PRN 26 from station CHAN.The vertical red arrow indicates the explosion time, 00:54 UT.The top right shows the enlarged TEC during 02:00-03:15 UT separately, to discern the arrival time of the TIDs.The capital letters "A," "B" and "C" mark the larger disturbances in STEC.Disturbance series extracted by the multi-step multi-order numerical difference method from station CHAN are presented in Fig. 4. As seen from Fig. 4a, TID "A" approximately 02:00 UT is discernable when the time step is set to 30 s, which is consistent with the results in Park et al. (2011).However, according to the top right of Fig. 3, the ionospheric STEC series also contains significant disturbances during 02:30-03:15 UT (namely, TIDs "B" and "C"), suggesting that the 30 s step may miss important disturbance signals.After we extend the time step to 90 s, the two later TID events are detectable, while TID "A" is "missing" (Fig. 4b).No more significant TID events are detected if we increase the time step to 300 s (Fig. 4c).
To obtain the periods of the TIDs, the corresponding results of spectral analysis for disturbance series are shown in Fig. 4d-f.Figure 4d shows that TID "A" has a frequency spectrum of 11-15 mHz, corresponding to the period of 1.1-1.6 min.The TIDs "B" and "C" have similar frequency spectrums of 3.3-7.5 mHz (corresponding to periods of 2.2-5 min) in panel e.The TID signal in Fig. 4f is not evident and the frequency spectrum of 1 to 2 mHz belongs to the ambient noise signal, which is the residual term when the ionospheric TEC is detrended.Based on the above results, two types of TIDs are observed by the GPS station CHAN after the North Korean UNE: TID "A" is a high-frequency wave with periods of approximately 1 to 2 min, whereas TIDs "B" and "C" are low-frequency waves with periods of approximately 2 to 5 min.
The two types of TIDs are also observed by other stations.Figures 5 and 6 show all the extracted disturbance series containing high-frequency TIDs and low-frequency TIDs, respectively; and corresponding attributes of the TIDs are listed in Tables 1 and 2, respectively.As a comparison, 3-day data centered on 25 May 2009 are illustrated both in Figs. 5 and  6.According to the results, the disturbance signals are significant on event day with larger signal to noise rate (SNR) varying from 5.4 to 33.6, but not on other two days.The disturbance signal is like an N wave with short duration, which is also observed on the ionospheric TEC after earthquake (Afraimovich et al., 2001).As shown in Tables 1 and  2, the periods for high-frequency TIDs are approximately 1 to 2 min and the low-frequency TIDs have periods about 2-5 min.Small differences in detected periods exist between different stations, due to the Doppler effect caused by the relative motion between the GPS satellite and the TIDs (Wang et al., 2007;Garrison et al., 2007).Both the frequencies of high-frequency and low-frequency TIDs are larger than the acoustic cutoff frequency (about 3.3 mHz in the lower atmosphere; Artu et al., 2004), suggesting that the two types of TIDs can be attributed to the acoustic waves.The ionospheric pierce points (IPPs), where the detected TIDs have maximum amplitudes, are also listed in Tables 1 and 2 and are illustrated in Fig. 2. Obviously, both the high-frequency waves and low-frequency waves are centered on the epicenter of the UNE.As these TIDs surround the epicenter and have similar features, we confirm that they were induced by the UNE rather than by other random events.
The average propagation velocity of TIDs is obtained using the linear distance between epicenter and IPP divided by propagation time, which is approximately equal to the trace velocity.The results are listed in Tables 1 and 2 and are illustrated in Fig. 2 as well.According to Fig. 2, in general, TIDs near the epicenter have larger propagation velocities than TIDs away from the epicenter, suggesting the propagation velocities feature distance-dependent attenuation.This phenomenon is consistent with the view that these TIDs are related to the acoustic waves induced by the UNE.Generally, for such a scenario, the trace velocity is basically c • [sec(e) − 1], where c is the sound speed and e is the elevation of the ray path of a wave as it leaves the epicenter (Kinsler et al., 2000).In addition, the TIDs on the northern side of epicenter have significantly lower propagation velocities than the southern parts.The possible reason may be because the northward-propagating disturbances are sensitive to the effects of Earth's magnetic fields (Heki and Ping, 2005;Otsuka et al., 2006;Occhipinti et al., 2008;Liu et al., 2011).Another phenomenon should be noted that many observed TIDs have velocities below the lowest sound speed within Earth's atmosphere at these latitudes.This may be due to the effects of Earth's magnetic fields mentioned above.In addition, the real traveling distance of acoustic waves induced by UNE is generally longer than the distance between the IPP and epicenter (sometimes may be ducted within the thermosphere), so the velocity is underestimated.
Comparing the results between Figs. 5 and 6, the amplitude of high-frequency TIDs is generally lower than the lowfrequency TIDs: the former is about 0.1 TECU whereas the latter is generally above 0.2 TECU.A possible reason for this phenomenon is that the high-frequency signal only reaches the lower ionosphere due to the viscous dissipation within the thermosphere which is inversely proportional to frequency.The acoustic waves or acoustic-gravity waves in lower ionosphere can also induce detectable disturbances (Choosakul et al., 2009;Watada and Kanamori, 2010;Nina and Čadež, 2013).This explanation can be strengthened by the results: observed average velocities of high-frequency TIDs are basically larger than the low-frequency TIDs.As the real height is less than 350 km, the calculated average velocities of highfrequency TIDs are amplified.The most likely origin for the two different TID components is stratospheric filtering of the acoustic waves generated by the UNE.This wave spectrum can be filtered to some degree by the wind profile in the stratosphere, and is further altered by nonlinear interactions as the waves break in the upper atmosphere (Moo, 1976).
According to the observation results in Yang et al. (2012), TIDs with periods of 2 to 5 min were detected after the North Korean UNE in 2009, while the high-frequency TIDs with periods of 1 to 2 min were not reported in that paper.The most likely reason is that the high-frequency (1-2 min) signals are relatively weak and cannot be detected easily, which suggests our detection method is also effective for weaker signals.Furthermore, there are also very long period (3-12 min) disturbances in data from PRN 21 that were reported in Yang et al. (2012).However, the periods of the same disturbances were both 2.4 to 5.2 min when employing our method and the frequently used polynomial fitting method (Wang et al., 2007;Tsugawa et al., 2007;Galvan et al., 2011).In essence, our method and the polynomial fitting method obtain the TEC disturbances directly in the time domain, whereas the method of Yang et al. (2012) is in the frequency domain.Considering the complexity and uncertainty in selecting suitable filters and the range of acoustic cutoff frequency, our results may be more reliable.

Summaries
This study employs a multi-step multi-order numerical difference method to detect TIDs using GPS observations from several IGS stations located around the epicenter of the North Korean underground nuclear test on 25 May 2009.Compared to the current numerical difference method, the modified method has a wider period spectrum for TID detection and possesses the capacity to identify different types of ionospheric disturbances.The detection method is also sensitive to relatively weak signals due to its magnification effect.
Two types of ionospheric disturbances, namely, highfrequency and low-frequency TIDs, have been observed surrounding the epicenter.The high-frequency TIDs generally have periods of 1 to 2 min, while the period range for lowfrequency TIDs is between 2 to 5 min.The propagation velocities of TIDs exhibit distance-dependent attenuation and north-south asymmetry.The high-frequency and lowfrequency TIDs can be ascribed to the acoustic waves in lower ionosphere and higher ionosphere, which originated from the UNE event on 25 May 2009.In short, this work demonstrates the capability of using GPS observations to detect UNE-generated ionospheric disturbances.At present, thousands of GPS stations from IGS and other organizations have been deployed across the globe including local districts.These GPS observations can be used as an important data source to monitor nuclear weapons development.

Figure 1 .
Figure 1.Amplification factors of TIDs for different time steps.

Figure 2 .
Figure 2. The locations of GPS stations (black triangles) and IPPs of detected TIDs.The black star represents the epicenter location.The circles denote the high-frequency TIDs and the squares denote the low-frequency TIDs.The color bar indicates the TIDs' propagation velocities.

Figure 3 .
Figure 3. Ionospheric STEC obtained by GPS observations of satellite PRN 26 at station CHAN.The three arrows in top right indicate the arrival time of TIDs.

Figure 4 .
Figure 4.The disturbance series extracted by multi-step multi-order numerical difference method at station CHAN.The order is set to 2. Panels (a-c) present the results of time steps 30, 90 and 300 s, respectively.Panels (d-f) are the corresponding spectrograms.

Figure 5 .
Figure5.The extracted disturbance series containing the highfrequency TIDs at different stations for 3-day data.The abscissa represents the UT time and the ordinate represents the disturbance series with unit of TECU.The title for every panel is composed of station ID and satellite PRN number.The bold red line indicates the series on event day and the green and blue for the former day and the latter day, respectively.

Figure 6 .
Figure6.The extracted disturbance series containing the lowfrequency TIDs at different stations for 3-day data.The abscissa represents the UT time and the ordinate represents the disturbance series with unit of TECU.The title for every panel is composed of station ID and satellite PRN number.The bold red line indicates the series on event day and the green and blue for the former day and the latter day, respectively.

Table 1 .
The propagation characteristics of observed highfrequency TIDs in terms of location (latitude and longitude), period and average propagation velocity.The SNR is the signal and noise rate.

Table 2 .
Same as Table 1 but for the low-frequency TIDs.