Remote Sensing Small Explosives With an Ionospheric Radar

Earth's ionosphere has long been targeted as a medium for remote sensing of explosive terrestrial events such as earthquakes, volcanic eruptions, and nuclear/conventional weapon detonations. Until now, the only confirmed ionospheric detections have been of very large events that were easily detectable through other traditional global sensor systems (e.g., seismic). We present the first clear, confirmed detections of relatively low yield 1‐ton TNT‐equivalent chemical explosions using pulsed Doppler radar observations of isodensity layers in the ionospheric E region. The shape and spectra of the detected waveforms closely match predictions from the acoustic ray tracing and weakly nonlinear waveform propagation models. The explosions were roughly three orders of magnitude lower yield than any previous confirmed ionospheric detection and represent the first conclusive evidence that explosions of this size can have clear impacts on the ionosphere. This technique could improve the remote detection of both anthropogenic and natural explosive events.


Introduction
Atmospheric waves, produced by energetic events at the Earth's surface, underground, or in the atmosphere, can propagate up to the thermosphere/ionosphere and couple with ionospheric plasma (Astafyeva, 2019;Huang et al., 2019;Meng et al., 2019). Propagation to these heights is advantageous as the free electrons in the ionosphere are an easy target for radio frequency (RF) remote sensing. While typical infrasound sensors such as microbarometers require direct contact with the wave, an RF sensor on the ground has a view of the ionosphere out to a few thousand km radius. RF remote sensing has the potential to detect events out to great distances, where one only needs to wait until the infrasound reaches the ionosphere rather than directly contacting the sensor. This technique to identify signals of interest has been successfully applied to volcanic eruptions (Themens et al., 2022), earthquakes (Maruyama et al., 2016), thunderstorms (Lay et al., 2015), rocket launches (Mabie et al., 2016), and large (∼kiloton TNT equivalent) nuclear and conventional explosives (Barry et al., 1966;Blanc & Rickel, 1989;Jacobson et al., 1988;Kundu et al., 2021;Park et al., 2011).
Most detections of ionospheric fluctuations caused by terrestrial events are carried out using global navigation satellite system (GNSS) receivers, where impulsive events cause detectable fluctuations in total electron content (TEC). For instance, the disastrous 4 August 2020 ∼1 kT explosion in Beirut was detected using such a method (Kundu et al., 2021). However, this detection was hardly obvious in the data, as other fluctuations of similar amplitude occurred well before the explosion occurred. The event did not appear to make a significantly large fluctuation in TEC despite being one of the largest chemical explosions to ever occur. This is likely due to the fact that the wavelengths of the infrasound generated by the explosion were significantly smaller than the integrated lengths of the ionosphere, meaning that much of the signal would be integrated out.
Techniques such as medium frequency (MF, 0.3-3 MHz) and high frequency (HF, 3-30 MHz) Doppler radar probe a single isodensity layer in the ionosphere, and therefore are not subject to signal loss through integration. Barry et al. (1966), successfully measured the signature of a 500-ton explosion in Alberta CA using a pulsed HF Doppler radar. About 9 min after the explosion, the shock wave clearly perturbed an F layer isodensity by 100 radians with a period of about 2 min. Similarly Blanc and Rickel (1989) describes an experiment where a 5 kt chemical explosion was detected in the ionosphere using a vertical incidence ionosonde and Doppler radar. The authors noted visible changes in ionospheric heights and the E region Doppler spectrum that were possibly caused by the passing of the shock wave.
While it is technically infeasible for the authors of this paper to reproduce a ∼kT chemical explosion today, ionospheric detections of much lower yield events are still highly relevant. It is of particular interest to determine the range of events that can and cannot be detected. Certainly the total energy is an important factor for detection. Higher energy events are more likely to have a larger impact. Perhaps equally important, however, are the temporal and spatial scales of the events, which dictate the spectral content. For instance higher frequencies are subjected to more attenuation, meaning that lower frequencies can propagate to higher altitudes. Additionally, as the wavefront propagates, nonlinear period lengthening occurs, where the high pressure travels faster than the low pressure. Therefore, the peak wavelength (frequency) measured increases (decreases) with increasing altitude. Fitzgerald and Carlos (1997) made continuous wave (CW) Doppler measurements of the ionospheric E region during a period when numerous half-ton TNT equivalent tests were being conducted nearby. The shots do not appear have been coordinated, and the exact timing and location of the shots was not known. However, during the time period in question, the authors found numerous Doppler chirps lasting a few seconds each. The authors interpreted the chirps as being signatures of explosions triggered by momentary instances of diffraction caused by laterally propagating perturbations to the reflecting isodensity surface. It is possible that these low yield detonations did indeed cause the detected signals, but without exact timing and location information it is difficult to say for certain.
Recently, in Obenberger et al. (2022) we presented two experiments to detect small 1-ton TNT equivalent explosions using Doppler sounding of the lower ionospheric E region. Shots of this magnitude have spectral content peaking near a few Hz, with no measured content below 0.1 Hz . Modeling the waveform using the InfraGA/GeoAc package provided by Los Alamos National Laboratory, we found that the signal would have a peak frequency between 0.9 and 0.4 Hz in at the altitude range between 96 and 110 km. Moreover, we found that the amplitude would decrease significantly over this range, implying that lower altitudes would be the best to target. For more information on InfraGA see Blom (2019) and Blom and Waxler (2012, 2017. Of the two experiments presented in Obenberger et al. (2022), only one resulted in detections, and these were low SNR and somewhat perplexing. The arrival time and spectral content fit well to the modeled expectation, but the temporal characteristics were difficult to interpret. Each shot had a different appearance and signal strength despite being of identical yield and temporally separated by 90 s.
In Obenberger et al. (2022), the closest ionospheric measurement point was 70 km horizontal from the test site. Infrasound modeling, however, consistently shows that the strongest thermospheric signals should be directly above the explosion site (e.g., Sabatini et al., 2019). We therefore have recreated our previous experiments using one transmitter and three receivers each with measurement points (mid-point of bistatic radar) within 25 km horizontal distance from the explosion site. We present the results from this experiment in this paper.

Experimental Setup
Humming Lobo was an experiment carried out by Sandia National Laboratories at the Energetic Materials Research and Testing Center (EMRTC) in Socorro NM on 28 March 2022. The campaign's primary scientific goal is to evaluate changes in the stratospheric infrasound propagation regime across variable seasons. Detonations were repeated at 180 s (3 min) and several hour intervals to examine minute and hour-scale variability in seismoacoustic signal propagation. The experiment was recorded across a network of local distance (<25 km) acoustic sensors, regional distance (100-400 km) co-located seismic and acoustic sensors and airborne balloon-based infrasound sensors. The experiment consisted of two sets of two 1-ton TNT equivalent chemical explosions. In each set the two shots were separated by 3 min.
The Air Force Research Laboratory (AFRL) participated in Humming Lobo by deploying a Digisonde Portable Sounder 4D (DPS4D) at White Sands Missile Range along with three MF/HF receive sites near Socorro NM. The geometry was set up such that the mid-points between the DPS4D and each of the receivers formed a triangle around EMRTC. The DPS4D was run at two frequencies separated by 200 kHz, namely 2.4 and 2.6 MHz in O mode only (X experiences more absorption). This mode allowed us to monitor the vertical movements of ionospheric layers at ranges with these particular plasma frequencies. For more information of the DPS4D see Reinisch et al. (2008Reinisch et al. ( , 2009. Figure 1 shows a map and lists the coordinates of the EMRTC test site, DPS4D transmitter (TX), MF/HF receiver (RX) sites and their respective mid-points. We also included a 25 km radius circle to show that the RX/TX links are all probing ionospheric points within 25 km horizontal distance from the test site. For this experiment the northernmost receive site was located at the Long Wavelength Array (LWA) station at Sevilleta National Wildlife Refuge (LWA-SV). The westernmost was located at the LWA station near the North Arm of the Very Large Array (LWA-NA) and the southernmost was located at a bunker at EMRTC just over 2 km from the test site. Figure 1 also shows a schematic of the vertical RF Doppler signal path and the movement of the reflective isodensity. As the infrasound from the explosion reaches the otherwise stationary isodensity surface, it displaces the surface from its resting position. This vertical displacement is detected as an oscillating Doppler signal. This creates what is effectively a RF microphone, where the isodensity surface is the diaphragm and the radar detects its vibrations. In reality the isodensity surface is also perturbed by numerous natural drivers creating a noise background for the radar.
Each RX used in this experiment is comprised of a simple cross-loop antenna connected to a two channel Ettus X300 universal software radio peripheral (USRP) with the basic RX daughter card, and Linux PC with a large storage device for recording data. Additionally, before each USRP there was a 30 dB low noise amplifier, and a 25-30 dB pre-amplifier at the antenna. At LWA-SV we used 30 dB pre-amplifiers and the other two RX , medium frequency/high frequency receiver (RX) sites, their respective mid-points, and a red circle with a 25 km radius from the test site. Also, inset on the bottom left corner is a schematic showing the vertical geometry of the Doppler signal path (TX to RX), the infrasound propagating from the explosion, the ionospheric isodensity, and its displacement caused by the infrasound. We note that the drawing is not done to scale.
10.1029/2022EA002737 4 of 8 locations used 25 dB units. For more information on the RX processing see Obenberger et al. (2022). Once range gated, the primary observable for each RX is the phase (ϕ) of the signal at each frequency. The phase is converted to Doppler through the time derivative (dϕ/dt). As a plasma isodensity in the ionosphere is perturbed by an infrasound wave the signature will be captured in the Doppler signal. We note that at the LWA-SV site we were actually running two separate X300 receivers, the antennas of which were separated by roughly 50 m.
The first shot in the first set occurred at 17:33:00.3 UTC and the second shot at 17:36:00.3 UTC on 28 March 2022. During this time there was good ionospheric propagation at the transmit frequencies and recent ionograms combined with the group delay indicated the reflection height to be near 100 km altitude. Furthermore, unlike the experiments presented in Obenberger et al. (2022), there was very little scintillation impact within an hour window of the first set. During the second set at 19:18:00.0 UTC, however, there was significant RF absorption due to a C9.8 Solar X-ray flare and subsequent enhanced D region (GOES-R Series Program, 2019). The signal strength during this time was too low to get accurate an Doppler measurement with the DPS4D so we will only present data from the first set.
Similar to  we can estimate the arrival times for each shot to the mid-points using InfraGA (Blom, 2019;Blom & Waxler, 2012, 2017. For an altitude of 100 km the mid-points for the LWA-SV, EMRTC, and LWA-NA RXs have arrival times of 5.52, 5.57, and 5.58 min, respectively. These estimates use temperatures and neutral densities from the Mass Spectrometer and Incoherent Scatter Radar Exosphere (MSISE) model and horizontal winds derived from the Horizontal Wind Model (HWM). For a description of HWM see Drob et al. (2015) and for MSISE see Picone et al. (2002).
We can also estimate the waveforms using the weakly nonlinear waveform propagator within the InfraGA model. Since the shots were of the same yield (1-ton TNT equivalent) used in Obenberger et al. (2022) we can assume that the input waveform used in that study is still valid. We do not know the exact altitude of the reflection points, but can infer from ionograms that they are in the lower E region between 95 and 110 km. Using this range we propagated the input waveform to these altitudes for the LWA-SV RX-TX mid-point. We can then use a continuous wavelet transform (CWT) of the pressure time series to analyze the peak frequency and aid visualization. As an example, Figure 2 shows the InfraGA modeled output using two shots separated by 3 min propagated to 100 km. Gaussian noise was added to simulate real data. The signals have a delay of 5.52 min and a peak frequency of 0.6 Hz. Similarly, for 95-110 km we find a range of 0.7 and 0.4 Hz, respectively. For the CWT analysis presented in this study we used a Morse wavelet and maintained the default MATLAB symmetry and time-bandwidth parameters, but used a upper frequency limit of 10 Hz and lower limit of 0 Hz. While the model presented here represents the expected pressure fluctuations of the neutrals at 100 km, the plasma is strongly coupled to the neutrals in the E region and should follow a similar temporal and spectral appearance to the neutral pressure wave. We note however, that the Doppler radar is really only measuring the phase changes along the propagation path to the reflecting layer, and these changes may be caused by additional effects other than simple displacement of the reflective isodensity. For instance Fitzgerald and Carlos (1997) suggest diffraction effects caused by lateral deviations in the reflective surface give rise to a Doppler signal. This is reasonable since the wavelength of the infrasound pulse is similar to that of the radar.

Results
Each RX measured the phase and amplitude (and therefore the Doppler) of the 2.4 and 2.6 MHz propagation channels at 50 Hz. To aid in the visualization of the Doppler signatures we use a CWT, which highlights spectral features as they evolve with time. Figure 3 shows the averaged CWT of the Doppler signal measured from the two receivers at LWA-SV and the CWTs from the single receivers at EMRTC and LWA-NA at 2.6 MHz. The CWT plots are zoomed in to show from 4 to 10 min after the first shot. It can clearly be seen that at ∼0.5 Hz two pulses were detected near 5.5 and 8.5 min at each RX location. While the properties of the pulses varied between each of the three stations, the two shots as measured by each respective station were quite consistent. At LWA-SV the peak frequency was near 0.5 Hz with an amplitude of 7.5 ± 0.6 m/s. EMRTC measured a slightly higher peak frequency at 0.71 Hz and a much higher amplitude of 24 ± 1 m/s. LWA-NA measured a frequency of 0.7 Hz with a peak amplitude of 7.5 ± 1 m/s. Despite some differences between each RX measurement, the fact that the arrival times closely match the InfraGA modeled expectations and the 3 min separation between pulses clearly indicates that the pulses were caused by the explosions. Indeed, the CWTs (especially from LWA-SV RX) look nearly identical to the InfraGA modeled expectations shown in Figure 2.
Figure 3 also shows the displacement (in meters) of the DPS4D signal as measured by each of the three stations near about 5.5 min after each shot. The displacement (d) is calculated simply by running the unwrapped phase (ϕ) through a high pass filter to remove larger scale trends and then using the simple relation d = λϕ/2π. Where λ is the wavelength of the radar. We note that the noise is generally higher in the LWA-NA and EMRTC measurements. The primary reason for this is that while the LWA-SV site benefited from two receivers whose CWTs were averaged together, LWA-NA and the EMRTC sites only had a single receiver each.
The SNR at 2.6 MHz was significantly higher than at 2.4 MHz for all three RX sites. While there would be slightly more ionospheric absorption at 2.4 MHz, lower RX and TX efficiencies at 2.4 MHz is likely to blame. For this reason we have primarily focused on 2.6 MHz. However, included as Supporting Information S1 are the CWTs for 2.4 MHz at each of the three RX sites. As can be seen when comparing to Figure 3, each has significantly lower SNR than at 2.6 MHz, and there are several scintillation spikes adding to confusion in the plots. As described in Obenberger et al. (2022), when a scintillation dip lowers the signal below the background noise level, the phase becomes completely random and therefore shows up as a broadband spike in the CWT. Despite this, we can see that the arrival times are approximately 3 s earlier than at 2.6 MHz, indicating a ∼1 km lower altitude for the 2.4 MHz isodensity.
Additionally, provided in Supporting Information S1 are CWTs of the amplitude time series from each RX site at both 2.6 and 2.4 MHz. As can be seen the pulses are easily identifiable in this manner as well. This implies that some amount of diffraction or interference is caused by interaction of the perturbed ionosphere and the RF signal. As the wavefront propagates laterally along the isodensity a changing diffraction pattern should arise as described in Fitzgerald and Carlos (1997).

Discussion
As the InfraGA model suggested, the measured arrival times are slightly different for each RX location. The pulses from each shot appear at LWA-SV, EMRTC, and LWA-NA near 5.35, 5.43, and 5.53 min, respectively. The arrival times are measured from the apparent onset of the displacement waveform as measured in Figure 3. Not surprisingly there is some disagreement with the InfraGA model, which is only as good as the input atmosphere. Differences between reality and the HWM and MSISE modeled atmosphere are likely to fault. Moreover, while the pulses arrived at nearly the same time after each shot as measured by the LWA-NA site, there is some discrepancy between the arrival times of the two shots as measured by the LWA-SV and EMRTC sites. At LWA-SV the second shot peaked 0.4 s later than the first shot, whereas at EMRTC the second shot peaked 0.7 s earlier than the first shot. Changing wind profiles, passing gravity waves, and a dynamic ionosphere all likely contributed to the sub-second difference in time of arrival measured for each shot. In Supporting Information S1 we have included plots of the temperature and wind profiles along with measurements made from radiosondes launched in Albuquerque NM at 00:00 and 12:00 UT. The radiosondes are able capture the temperatures and winds up to about 35 km altitude. We also include 3-hr assimilated meteorological estimates taken from the Modern-Era Retrospective analysis for Research and Applications version 2 (MERRA-2) from 15:00 UT. This data driven model provides parameters up to ∼65 km. These data agree relatively well to MSISE and HWM up to the maximum radiosonde height, above this more significant deviations occur between HWM and the MERRA-2 data.
The waveform shapes measured by each RX site are self-consistent between the two shots. However, it is intriguing that despite the fact that all three sites were probing ionospheric regions at about the same horizontal distance from the test site, we find some stark differences in waveform appearance. Here we can use the displacements shown in Figure 3 to compare the movement of the isodensity at each measurement point. As can be seen, the waveforms measured by each of the three sites show key differences. Perhaps the oddest in appearance is that measured by the EMRTC RX which is asymmetric and displays a quick positive to negative jump. Moreover, the EMRTC RX measured velocities that were more than double those measured by the other two sites. It is not immediately clear why the waveforms would be so different from location to location, but the mountainous terrain where the shots occurred along with a dynamic atmosphere likely played a role. Both the terrain and the perturbations on the isodensity are all of size scales that would require full wave modeling, and the ray tracing methodology used in InfraGA is not be able to capture such effects. The differences could also be the result of the radar interactions with the perturbed surface. The size of these perturbations are on the order of the maximum Fresnel zone of the radar, which means that full wave RF modeling is needed to fully understand the Doppler radar measurements.
Our results are consistent with the conclusions of Fitzgerald and Carlos (1997), as the events reported there were likely of similar yield and therefore should have been detected by a Doppler radar. However, without timing or location information, it is difficult to say for certain what exactly Fitzgerald and Carlos (1997) detected. Furthermore, our results are in agreement with (Barry et al., 1966); it appears that the Doppler radar technique is more sensitive than GNSS TEC based on comparisons to Kundu et al. (2021) and Park et al. (2011), which show only weak detections of massive kiloton explosions. As mentioned earlier, TEC measurements are inherently weaker because they integrate through most of the infrasound waveform. Additionally, we note that due to attenuation and non-linear period lengthening the amplitude of the infrasound decreases dramatically with increasing altitude. Simulations with InfraGA consistently show a rapidly decreasing amplitude of infrasound as it travels to higher altitudes.
The results presented in this paper indicate that explosives as small as 1-ton TNT equivalent produce measurable disturbances in the ionosphere. Explosions of this yield are relatively common in mining operations which can utilize yields ranging from grams up to even kilotons. While common in mining operations, such explosions can be devastating when used in terrorist attacks or war-zones. For scale the truck bomb used in the 1995 Oklahoma City bombing was estimated to be roughly 2.3 tons. Furthermore, a typical cruise missile can carry up to a half ton of conventional explosives, and aircraft dropped ordinance can range from 10s of kg to over 10 tons. It is safe to assume that in active war-zones ionospheric effects could be measured with an MF/HF radar, and a disturbed ionosphere may be possible. Further modeling is required to fully understand such effects.
In addition to enhanced remote sensing of small anthropogenic explosives, the ability to record the ionospheric imprint of much smaller events has favorable implications for remote sensing of geophysical phenomena on both Earth and other planets. The utility of infrasound recording using ground pressure sensors on Earth is well established, but this modality suffers from high latency due to the relatively slow lateral wave propagation speeds. In addition, wind noise can obscure faint acoustic arrivals from distant sources, and atmospheric structure can prevent the transmission of sound waves from the source to the receiver. This is a particular concern in areas like the Aleutian Islands, where distant arrays of infrasound sensors listen for explosions that may impact commercial flight lanes (Iezzi et al., 2019;Lamb et al., 2015). Continuous Doppler sounding of the ionosphere above targets of interest provides a lower latency monitoring method with less dependency on lower atmospheric conditions. Some planetary bodies are very difficult to monitor acoustically, either because they have an inhospitable surface (e.g., Venus) or lack a solid interface entirely (e.g., Jupiter). Ionospheric Doppler sounders orbiting Venus could provide information on geophysical activity of interest, such as volcanic eruptions and earthquakes, with better resolution than previously proposed airglow methods (Byrne & Krishnamoorthy, 2022).
Further development of this technique is critical before it could be used reliably in an operational sense. Future experiments would be better off utilizing a much more powerful TX source. This experiment utilized a pulsed DPS4D system with 300 W transmit power with a 5% duty cycle. It would be advantageous to use a higher power transmitter with frequency modulated continuous wave capability. Powerful transmitters such as coastal radars or AM radio stations could easily be used as signals of opportunity. It would also be enlightening to use a range of frequencies in order to compare the signal strength as a function of altitude. Additionally, future experiments should test a range of shot yields, across multiple seasons, and utilize more oblique propagation angles to investigate the reliability of this technique at great horizontal distances. While the demonstration presented in this paper utilizes Doppler sounding of near vertical incidence, our technique could easily be adapted for oblique cases targeting regions hundreds of km away. The only requirement would be that a TX to RX mid-point could be established nearly vertical to the target of interest. Or in the case of a powerful over the horizon radar, Doppler variations of ground scatter could be used for detection over a large area. Finally, both this study and Obenberger et al. (2022) indicate that even explosions taking place close in space and time are not always detected consistently. A comprehensive understanding of the prevalence and cause of such false negatives is urgently needed before this technique can be used in operational monitoring systems.

Conclusion
In this paper we have presented clear detections of ionospheric impacts from a set of two 1-ton TNT equivalent surface explosions. This study proves that the ionosphere can be used to detect explosions of this modest yield. The technique is simple, and utilizes an MF/HF Doppler radar to probe an isodensity in the lower E region. Essentially we form a RF microphone, where the isodensity surface acts like the diaphragm and the Doppler radar detects its vibrations. The detections presented in this paper represent the lowest yield explosions confirmed to be detected through ionospheric remote sensing and are nearly three orders of magnitude lower than previous confirmed successes.

Data Availability Statement
All of the data used in his paper can be found at the following links. The GOES data used in this paper can be found at https://www.ngdc.noaa.gov/stp/satellite/goes-r.html. The radar processed MF/HF Doppler data used to make the measurements reported in this paper are available at https://doi.org/10.5281/zenodo.7232404.