Atmospheric wave energy of the 2020 August 4 explosion in Beirut, Lebanon, from ionospheric disturbances

Atmospheric waves excited by strong surface explosions, both natural and anthropogenic, often disturb upper atmosphere. In this letter, we report an N-shaped pulse with period ~ 1.3 min propagating southward at ~ 0.8 km/s, observed as changes in ionospheric total electron content using continuous GNSS stations in Israel and Palestine, ~ 10 min after the August 4, 2020 chemical explosion in Beirut, Lebanon. The peak-to-peak amplitude of the disturbance reached ~ 2% of the background electrons, comparable to recently recorded volcanic explosions in the Japanese Islands. We also succeeded in reproducing the observed disturbances assuming acoustic waves propagating upward and their interaction with geomagnetic fields.

The Beirut explosion occurred around the sunset time, when strong ionospheric irregularities known as equatorial plasma bubbles (EPBs) often develop due to the Rayleigh-Taylor plasma instability and mask subtle changes in ionosphere 14,15 . Gentile et al. 16 showed that rates of EPBs occurrence depends on seasons and regions ( Supplementary Fig. 1), i.e. the high rate concentrates in spring and autumn worldwide and also in winter in the America-Atlantic-Africa region. The EPB production rate is not high at the longitude of Beirut in early August. We also note that geomagnetic activity was low around the time of the Beirut explosion as suggested by various indices (e.g. Kp = 1.7) (Supplementary Fig. 3). We did not find EPB signatures in the VTEC around the explosion time ( Supplementary Fig. 3), although moderate ionospheric scintillation signals are seen a few hours after the explosion ( Supplementary Fig. 4).
We examined TEC data from 15 continuous GNSS stations in Israel/Palestine. Figure 2a represents vertical TEC time series obtained using the GPS satellite 22 at five ground GNSS stations to the south of Beirut. Figure 2b shows the trajectory of sub-ionospheric points (SIPs) calculated assuming the thin ionosphere at altitude 300 km. Clear N-shaped disturbances emerge 10-12 min after the explosion. The signals are invisible at stations too close to Beirut (e.g. stations to the north of drag in Fig. 2). The signals appear to arrive later with smaller amplitudes at stations farther from Beirut. We could not find clear signals for other GPS satellites ( Supplementary Fig. 3). These GNSS stations also tracked GLONASS satellites, but we could not find clear signals because of unfavourable satellite distribution at the explosion time. We also note that such clear disturbance signals are either missing or ambiguous at GNSS stations far to the north, e.g. in Turkey and Cyprus. As shown in Heki 9 , such N-shaped signatures in TEC appear on the southern side of the explosion site in the northern hemisphere. Although the waves in neutral atmosphere propagates without such directivity, the electron movements are constrained along the geomagnetic field causing such directivity 17 . The disturbance amplitudes are also influenced by the difference in the angle between the wavefront and line-of-sight (LOS) of the satellites. These points are further discussed in the next section.
In Fig. 3, we drew a diagram indicating the distance between sub-ionospheric points (SIP) and Beirut as a function of time with colours showing the vertical TEC anomaly. There we used all the available GNSS station located to the south of Beirut. We found that the anomaly propagated southward at an apparent speed ~ 0.8 km/s. This is the sound wave velocity in the lower part of the F region of the ionosphere and much slower than the Rayleigh wave (~ 3.8 km/s) often found for coseismic ionospheric disturbances 18 , and much faster than the internal gravity wave (0.2-0.3 km/s) often excited by very large earthquakes 19 . Positive TEC anomalies seen around the explosion time within ~ 130 km from Beirut (Fig. 3) do not propagate from Beirut and would not be related to the explosion.
Comparison with numerical simulations. The acoustic waves propagate upward being refracted following the velocity structure shown in Fig. 4a, and waves emitted in high angles (within ~ 20° from zenith) reach and disturb the ionospheric F-region. Here, we try to reproduce the arrival times, waveforms, and relative intensities of the observed ionospheric disturbances made by the Beirut explosion with a simple numerical Sequence of the ash cloud captured during the explosion (snaps archived at https ://www.thene wsmin ute.com/artic le/shock ing-video s-show-power ful-explo sion-rocks -leban on-capit al-beiru t-13004 1). SkySat imagery before and after the explosion shows the impact of explosion in Beirut (imagery captured on May 31, 2020 and on August 5, 2020) (archived at the largest earth-observationsatellite-network of global dataset platform at https ://www.thegu ardia n.com/world /2020/aug/06/beiru t-explo sion-befor e-and-after -satel lite-image s). (b) Illustration of ionospheric disturbance caused by an explosion can be detected by differential ionospheric delays of microwave signals of two carrier frequencies from GNSS satellites. The disturbances in ionosphere are observed 10-12 min after the explosion, the time necessary for the acoustic wave to reach the ionospheric F region. The wave makes electron density anomalies (a pair of positive and negative anomalies shown with red and blue) on the southern side of the explosion. Figure (a www.nature.com/scientificreports/ simulation. Figure 4b shows the north-south vertical cross section of a space where an N-shaped acoustic pulse propagates upward. The period of the pulse was ~ 80 s (see Discussion), the period often observed in ionospheric disturbances by volcanic explosions 7,10 . Figure 4b also shows the position of the N-shaped acoustic pulse at three epochs, together with the LOSs connecting the GPS satellite 22 and five GNSS stations with clear explosion signatures in TEC. These LOSs first encounter the positive electron density anomaly 10-12 min after explosion, marking the onset of the TEC perturbation. Then LOS penetrate the negative anomaly making the TEC drops following the positive peak.
In illustrating the electron density anomalies made by the acoustic pulse, we assumed two attenuating factors, (a) height-dependence of the background electron density and (b) angle between the particle motion of neutral atmosphere and local geomagnetic field. For (a), we suppressed the electron density anomalies using a factor indicating the ratio of the electron density at that height relative to its maximum. We assumed the Chapman distribution 20 , of the altitude (z) dependence of the electron density,N(z) = N c exp 1 − ϕ − exp −ϕ /2 . There, ϕ = (z = h c )/h and the altitude of largest electron density h c was set to 300 km, and h is assumed 65 km (Fig. 4a).  www.nature.com/scientificreports/ For (b), we attenuated the anomalies by multiplying with directional cosines of the wave propagation directions onto the geomagnetic field. Next, based on the LOS connecting the GPS satellite 22 and the five ground stations, we calculated the sum of the electron density anomaly at each intersection of the ray and LOS. By repeating the calculation every 30 s, we have synthesized slant TEC change time series and compared them with observations in Fig. 4c. The synthesized TEC are multiplied with an arbitrary factor to adjust the anomaly amplitudes (i.e. only relative amplitudes of the simulation are meaningful). Nevertheless, we can see the consistency in the ratios of TEC disturbance amplitudes among stations as well as arrival times and waveforms, between the synthesized and the observed disturbances. Although the source function has equal amounts of positive and negative anomalies, the disturbance signals are dominated by positive changes. This is because the negative electron density anomalies are partly cancelled by the penetration of the LOS through positive parts. Absence of signals at GNSS stations too close to Beirut would be due to the cancellation of positive and negative electron density anomalies caused by high-angle penetration of the LOSs with the wavefront.

Discussions
The period of the N-shaped disturbance (~ 1.3 min) corresponds to the shortest period of the atmospheric bandpass filter ( Supplementary Fig. 5), suggesting that the original energy of the atmospheric waves by the Beirut explosion concentrates in shorter-period components. The ionospheric disturbances caused by the explosion was 0.28 ± 0.03 TECU in peak-to-peak amplitude of slant TEC (calculated from the data of drag, klhv, yrcm). Supplementary Fig. 6 shows the fit of synthesized time series with the data from the two stations, drag and klhv (same as the top two in Fig. 4c). There, we changed the period of the source function and found 80 s (~ 1.3 min) best reproduces the observed waveform.
Ionospheric disturbances by five volcanic explosions of four active volcanoes, observed using GNSS-TEC data in Japan, were recently compiled 10 . The amplitudes of the slant TEC disturbance, normalized by background vertical TEC, would serve as a measure for the intensity of volcanic explosions. The background vertical TEC of the Beirut explosion is ~ 13.4 TECU according to the Global Ionospheric Map 21 (Supplementary Fig. 7), which www.nature.com/scientificreports/ makes the ratio ~ 2.1%. This is comparable to recent volcanic explosions in Japan studied in Cahyadi et al. 10 , and is slightly larger than the 2004/9/1 (20:02 local time) eruption of the Asama Volcano, Central Japan, whose ionospheric disturbances were found for the first time with GNSS by Heki 9 (Supplementary Fig. 8).
Ionospheric disturbances with peak-to-peak amplitude of ~ 0.03 TECU by a surface mine blast in Wyoming, USA in 1996 with 1.5 kt of ANFO explosives (similar in power to TNT) were detected with GPS receivers 8 . Despite a slightly larger amount of explosives, its TEC disturbance is only ~ 1/10 of the present case. Although the exact background VTEC at that time is unknown, the difference seems significant. We think this small amplitude partly comes from the larger incidence angle of the LOS connecting the ground stations in Wyoming and GPS satellite 6 with the wavefront (stations are located between the blast site and SIPs). It would also be partly due to the design of the mine blast done in a pit to fracture surface rocks, which is different from the Beirut explosion that occurred on an unguarded surface.

Materials and methods
GNSS and space weather data. To estimate TEC variation during August 4, 2020 Beirut explosion, we have collected RINEX GNSS data from stations in Israel and Palestine, operated and maintained by survey of Israel and archived at Scripps Orbit and Permanent Array Centre (SOPAC, http://sopac -old.ucsd.edu/dataB rowse r.shtml ). We use the standard sampling rate of 30 s in the daily GNSS file in order to quantify the TEC changes in the ionosphere during the Beirut explosion. We also studied data from GNSS stations in Turkey, Iraq, Cyprus, but did not find clear signals related to the explosion.
We used the Kp and Dst indices, F10.7 and changes in geomagnetic horizontal (H) and vertical (Z) components to evaluate geomagnetic activities during the studied day ( Supplementary Fig. 2). The first two indices are based on geomagnetic field variations, and F10.7 indicates the changes in solar extreme ultraviolet radiation. The data are archived from OMNIWeb (https ://omniw eb.gsfc.nasa.gov/form/dx1.html). The geomagnetic data are from the Tihany station, Hungary (THY, 46.9° N, 17.9° E) (https ://www.inter magne t.org/data-donne e/downl oad-eng.php).

GNSS-TEC processing strategy.
In this study, we used GPS satellites to capture ionospheric disturbances.
The phase difference between the two frequency, L 1 (~ 1.5 GHz) and L 2 (~ 1.2 GHz), of the microwave signals from satellites located ~ 20,000 km above the Earth's surface provide information on the ionospheric electrons integrated along the LOS called slant TEC (STEC). The methods in GNSS-TEC are described in detail in e.g., Calais et al. 8 and Heki 10 . We first remove ambiguities in carrier phase differences by letting them align with differential pseudo-ranges (codes). The observed STEC is a combination of the true TEC and satellite/receiver biases: We express TEC in TEC unit (TECU) (1TECU = 10 16 electrons/m 2 ), which is related to the delay as, where t is the difference in delay between the two frequencies.
We used satellite biases included in the header information of the Global Ionospheric Map files 21 and determined the receiver bias using the minimum scalloping technique 22 . STEC values are often converted to absolute vertical TEC (VTEC) values by removing the inter-frequency biases in GNSS receivers and satellites and dividing by the obliquity factor S(θ) 23 , that depends on the satellite elevation angle θ.
where the obliquity factor S(θ ) is defined as The parameter β is the incidence angle of the LOS with the ionosphere at altitude h, and we used the mean radius of the Earth (6,378 km) for R e .
Modeling of GNSS-TEC time series. The arrival times, waveforms, and relative intensities of the observed ionospheric disturbances by the Beirut explosion are reproduced by adjusting a simple numerical function following Heki 9 . The function is made of a set of positive and negative pulses (red curve in Fig. 4c) and expressed as where a represents amplitude of the disturbance. This numerical function has maximum and minimum at t = −σ and t = σ respectively, and we assumed that the explosion occurred at t = −2σ . We used σ = 20, which gives 80 s (~ 1.3 min) as the period of the disturbance. As shown in Supplementary Fig. 6, this value best reproduces the observed disturbances. This period approximately corresponds to the shortest period of the band-pass filter by the atmosphere (Supplementary Fig. 5).  www.nature.com/scientificreports/

Data availability
The GNSS and Space weather data used in this paper can be obtained online