Chemical explosions are ground truth events that provide data, which, in turn, can enhance the understanding of wave propagation, damage assessment, and yield estimation. On 4 August 2020, Beirut, Lebanon was shocked by a catastrophic explosion, which caused devastating damage to the Mediterranean city. A second strong chemical explosion took place at the Xiangshui, China chemical plant on 21 March 2019. Both events generated shock waves that transitioned to infrasound waves, seismic waves, as well as hydroacoustic signals with accompanying T-phases in the case of the Beirut event. In this work, the seismo-acoustic signatures, yields, and associated damage of the two events are investigated. The differentiainterferometry synthetic aperture radar analysis quantified the surface damage and the estimated yield range, equivalent to 2,4,6-trinitrotoluene [C7H5(NO2)3] (TNT), through a “boom” relation of the peak overpressure was evaluated. Infrasound propagation modeling identified a strong duct in the stratosphere with the propagation to the west in the case of the Beirut-Port explosion. In the case of the Xiangshui explosion, the modeling supports the tropospheric propagation toward the Kochi University of Technology (KUT) sensor network in Japan. Although the Beirut yield (202–270 ± 100 tons) was slightly larger than the Xiangshui yield (201 ± 83.5 tons), the near-source damage areas are almost the same based on the distribution of damaged buildings surrounding the explosions.

Industrial explosions can be tragic events, which leave devastation, loss, and casualties behind, but for forensic scientists, industrial explosions generate valuable data and records, which, when combined with ground truth information, can be used to mitigate the effects of future disasters. Seismo-acoustic and remote sensing data are valuable for calibrating propagation models, making yield estimates, and quantifying damage scaling. In this paper, we analyze the data from two large chemical explosions, which were recorded on seismic and infrasound sensors. Additionally, remote sensing data were used to quantify the explosive damage in the neighboring areas. This assessment can provide valuable information for monitoring artificial explosions because single-fired chemical explosions can be scaled and compared to the effects of nuclear explosions (Stump et al., 1999). The first event analyzed is the Beirut-Port (B-p) explosion, which occurred at 15:08:18 Greenwich Mean Time (GMT) on 4 August 2020. A large portion of the stored ammonium nitrate in the port detonated, leading to at least 200 deaths, 6500 injuries, an estimated 300 000 homeless, and capital loss of at least $15  × 109 (BBC, 2020). The explosion was widely documented and many witnesses recorded the fire, detonation, and consequent blast waves on video, and these were shared on social media (Hubbard, 2020; Livescience, 2020). Additionally, the explosion generated seismic, acoustic, infrasound, and hydroacoustic waves, which propagated through the lithosphere, atmosphere, and hydrosphere. The signals were widely recorded regionally.

The Xiangshui chemical plant explosion, which occurred in the afternoon of 21 March 2019, is the second event analyzed. A strong explosion took place in a fertilizer factory at Chenjiagang Chemical Industry Park, Chenjiagang, Xiangshui County, Yancheng, Jiangsu, China. At least 94 were severely injured and 78 people were killed and the factory was totally destroyed along with many neighboring facilities (BBC, 2019; Zhang et al., 2019). This explosion, as well, was recorded by seismic and infrasound stations. Open-source data from China were limited but strong infrasound was recorded in Japan. In this work, we document the seismo-acoustic signatures of the two events and estimate the sizes of the sources using different techniques.

A fire had started in hanger number 12 at the B-p, as evidenced by gray soot and smoke billowing from the hanger (Fig. 1). Ten minutes later, a bright burst, followed by sprays of smaller flashes, was seeen and attributed to being the fireworks stored in the hanger going off after catching fire. In the same hanger, 2.7 kt of ammonium nitrate had been stored for a long time (Hubbard, 2020). Ammonium nitrate is widely used as a fertilizer and in manufacturing explosives. Ammonium nitrate cannot ignite “on its own” and needs another explosion to initiate the reaction. It is believed that the detonation of the fireworks provided the condition for the transition from deflagration to detonation in the ammonium nitrate. Once initiated, the ammonium nitrate exploded violently, producing a strong flash as shown in Fig. 1(c). In the detonation process, solid ammonium nitrates rapidly decomposes into gaseous nitrous oxide and water vapor (AIGA, 2008). Water vapor can be clearly distinguished in Fig. 1(f) as the cloud color accompanying the shock wave with the nitrogen dioxides, forming the brownish gases that evaporate upward.

FIG. 1.

(Color online) The images depict (a) the beginning of a small fire a few seconds before the explosion, (b) the fire spreading, initiating the fireworks stored nearby, (c) the fire providing conditions for the explosion, and (d) the ammonium nitrate explosion. (e),(f) The generation of the fireball and initiation of the pressure waves is shown. (g) The shock waves reach the buildings near the explosion site. Panels are stills from Newsflare video “Extended slow-motion HD version of devastating Beirut blast” available at https://www.newsflare.com/video/372044/extended-slow-motion-hd-version-of-devastating-beirut-blast (Last viewed February 25, 2022). Credit: Agoston Nemeth / Newsflare.

FIG. 1.

(Color online) The images depict (a) the beginning of a small fire a few seconds before the explosion, (b) the fire spreading, initiating the fireworks stored nearby, (c) the fire providing conditions for the explosion, and (d) the ammonium nitrate explosion. (e),(f) The generation of the fireball and initiation of the pressure waves is shown. (g) The shock waves reach the buildings near the explosion site. Panels are stills from Newsflare video “Extended slow-motion HD version of devastating Beirut blast” available at https://www.newsflare.com/video/372044/extended-slow-motion-hd-version-of-devastating-beirut-blast (Last viewed February 25, 2022). Credit: Agoston Nemeth / Newsflare.

Close modal

The explosion signals were widely recorded by the seismic stations in Israel as well as those in the occupied land of Palestine, Jordan, Cyprus, and Egypt. The nearest seismic station was GEM, from the Israel National Seismic Network, at a distance of 77.7 km, and the farthest station with a clear seismic record was NDB3 of the Egyptian National Seismic Network at a distance of 727.1 km. The seismic records illustrate a collection of phases from the explosion with a P wave as the first arrival at all stations, followed by the hydroacoustic signals observed on the stations near the coastline and, last, acoustic signals. The hydroacoustic signals can be generated by the explosion and travel through the sea, whereas the acoustic signal may be refracted through the atmosphere. Table IV (see the  Appendix) lists the seismic stations and arrival times of the different phases from the recorded waveforms shown in Fig. 2.

FIG. 2.

(Color online) The seismic waveforms in the Eastern Mediterranean (Israel, Jordan, Cyprus, and Egypt) are shown.

FIG. 2.

(Color online) The seismic waveforms in the Eastern Mediterranean (Israel, Jordan, Cyprus, and Egypt) are shown.

Close modal

The hydroacoustic signals are recorded on coastal and ocean-bottom seismometers. The nearest station that recorded such a signal was CY606, an element of the Cyprian ocean-bottom seismic array, at a distance of 103.2 km, whereas the furthest station which recorded the signal was NDB3 at a distance of 727.1 km. The acoustic signal is recorded on each station with distances ranging from 212 to 316 km from the epicenter. There are no signals recorded on stations located at distances less than 200 km, which is the acoustic shadow zone. The recorded signals at distances farther than 316 km were too weak to be detected by the seismic sensors. The arrival times of the three phase types were poltted against the distance to estimate the phase velocity (Fig. 3). The arrivals were fit with a linear least square algorithm to estimate the average velocity of each phase. Figure 3 shows the fits and velocities of the different phases. The P wave velocity was 6.39 km/s, the hydroacoustic velocity was 1.51 km/s, and the acoustic velocity was 0.301 km/s.

FIG. 3.

(Color online) The plotted data and three fitting lines of the average velocities and different phase speeds are shown.

FIG. 3.

(Color online) The plotted data and three fitting lines of the average velocities and different phase speeds are shown.

Close modal

The seismic records can be used to calculate the yield of the explosion in terms of the equivalent measures of 2,4,6-trinitrotoluene [C7H5(NO2)3] (TNT) through the relation between the seismic magnitude and yield. The magnitude characterizes the relative size of the event based on the measurement of the maximum recorded motion. There are different types of magnitudes, including the local magnitude. For this event, we used the local magnitude (ML) to estimate the yield of the Beirut explosion. Hutton and Boore (1987) defined the local magnitude scale following the earlier work by Richter (1935), which uses the maximum trace amplitudes recorded on Wood-Anderson (WA) horizontal seismometers as

ML=logAmax(WA)logA0,
(1)
log(A0)=1.110logΔ100+0.00189(Δ100)+3.0,
(2)

where Amax(WA) is the maximum trace amplitudes (mm) recorded on a simulated WA horizontal seismometer [Eq. (1)]. logA0 is the station correction, which depends on the attenuation due to the path effect of the signal from the explosion site to the recorded seismic stations [Hutton and Boore, 1987; Eq. (2), where Δ is the distance in kilometers between the seismic site and the epicenter location]. The digital waveforms, from the 20 regional three-component seismic stations to distances of up to 527 km, were analyzed and processed to estimate the ML magnitude of the explosion using the Seismic Analysis Code (“SAC”; Goldstein and Snoke, 2005). Each signal is corrected for its instrumental response and then convoluted with WA's response (Fig. 6). Using the WA record, the maximum recorded displacement from each horizontal component is measured to estimate the average value. The original explosion signals and simulated WA seismograms are presented in Fig. 4 for the stations MMLI, Δ = 163 km [Fig. 4(a)]; SALP, Δ = 205 km [Fig. 4(b)]; and GHAJ, Δ = 289 km [Fig. 4(c)]. The estimated local magnitudes obtained from each of the 20 seismic stations range from 3.2 to 3.9 (see Table V in the  Appendix). The relation described by Archuleta et al. (1982) was used to determine the average values of the magnitudes from the 20 seismic stations,

xavg=10((1/N)i=1Nlogxi),
(3)

and the standard deviation,

SD[logxavg]=[1Ni=1N(logxilogxavg)2]1/2.
(4)
FIG. 4.

(Color online) The seismograms of the original signal recorded explosion and their simulated WA traces from SALP (a), GHAJ (b), and MMLI (c) seismic stations are shown.

FIG. 4.

(Color online) The seismograms of the original signal recorded explosion and their simulated WA traces from SALP (a), GHAJ (b), and MMLI (c) seismic stations are shown.

Close modal

The resulting local magnitude of the Lebanon explosion is 3.55 (± 0.15 ML).

There are many relations between the magnitude and equivalent charge weight (TNT equivalent) for different regions and different types of explosive sources. The most appropriate relations in this regard are the relations calculated for the Dead Sea region, which is near by. The relation is given by Gitterman et al. (2005) as: M=0.2937+0.7327log(W,kg). Application of this model produces an average yield estimate of 202.2 tons (± 127.55 tons).

Infrasound can be defined as an inaudible sound with a frequency range lower than 20 Hz (ElGabry et al., 2017). A large number of chemical and nuclear explosions occurred between 1970 and 1990 with yields up to 8 kt (1ton=103kg); however, the number of controlled experiments and explosions, in general, decreased in the last decades (Gitterman, 2010). The interaction of seismic and infrasound signals accompanying explosions can be complex, namely, seismic waves are generated from the direct coupling of the explosion to the ground and air-to-ground coupling. An explosion in or near water can generate hydroacoustic waves and, under some circumstances, the detectability of the explosions improves (Ceranna et al., 2009). Through the low-velocity sound channels, the acoustic signal can propagate to thousands of kilometers (Johnson, 2003). Stratospheric winds produce ducts for the propagation of low-frequency energy, which refracts back to Earth's surface. The arrivals, which refract from the stratosphere, can be detected at distances of over 200 km (Green et al., 2011). The effective sound speed for altitudes greater than 100 km increases in proportion to the temperature increase, suggesting mesospheric arrivals, although the amplitudes may be highly attenuated (Pilger et al., 2019). The simulation of the infrasound propagation in the atmosphere from the source to the receiver can be modeled by different methods: classical ray tracing (Georges, 1972), normal modes (Pierce, 1967), and the parabolic equations (PEs; Lingevitch et al., 2002). In the case of the assumption of the plane-parallel wavefront in a linear medium, classical ray tracing can be applied. Groves (1955) and Thompson (1972) enhanced the application of the classical ray tracing theory by neglecting the effect of vertical winds and assuming that the atmosphere is a horizontally stratified field (Drob et al., 2003). The PE technique is used to provide full wave, wind-dependent propagation estimates with high accuracy. With the application of the wide-angle assumption and high Mach number, the PE in the infrasonic propagation becomes more efficient (Lingevitch et al., 2002). In this study, the NCPAprop package was used for the infrasound propagation modeling (Waxler et al., 2017b). The B-p explosion was detected west of the explosion. The infrasound waves propagated for thousands of kilometers and were receded at I48TN, I26DE, and I17CI, in addition to the weak detections at the I42PT and I11CV stations. Furthermore, a single infrasound sensor of the Egyptian infrasound station HLW recorded a weak signal. The PMCC (progressive multichannel correlation) analyses of the recorded infrasound waveform data at I48TN, I26DE, and I17CI are used to estimate the azimuth and apparent velocity at each infrasound array as shown in Fig. 5. The PMCC was applied to 19 frequency bands from 0.1 to 7 Hz by identifying the frequency limits for each step with 1/3 octave bands and an overlap of 95% between windows lengths. In addition, we assumed the minimum apparent velocity of detection to be 300 m/s.

FIG. 5.

(Color online) The beam waveforms for the IMS stations (I48TN, I26D, and I17CI, respectively) using the PMCC are shown.

FIG. 5.

(Color online) The beam waveforms for the IMS stations (I48TN, I26D, and I17CI, respectively) using the PMCC are shown.

Close modal

The apparent velocities for these sites vary from 340 to 360 m/s as shown in Table I. These velocities reflect the effects of the wind and propagation path. The explosion location is reflected by the calculated back-azimuths at the different sites. The periodic patterns observed along the azimuth to I48TN [Fig. 9(a)] are consistent with multiple reflections as predicted in the propagation modeling in Fig. 8. The observed celerities from the stations are well-correlated with the estimated values from the ray tracing modeling as shown in Table I. The celerities computed from the upper stratosphere range from 295 to 310 m/s.

TABLE I.

The infrasonic arrivals to the International Monitoring System (IMS) infrasound stations.

DurationBack-azimuthApparent velocity rangeFrequency rangeObserved celerityComputed celerity
StationTime(min)(deg)(m/s)(Hz)(m/s)(m/s)
I48TN 17:06:41 21:13.1 81.5–93.0 336.0–392.0 0.2–6.4 285.3–336.4 297 
I26DE 17:12:21 15:25.8 119.2–134.4 328.0–401.0 0.1–2.0 292.7–329.1 297 
I17CI 19:44:32 13:33.1 44.0–50.0 302.0–364.0 0.1–1.6 293.6–308.1 314 
DurationBack-azimuthApparent velocity rangeFrequency rangeObserved celerityComputed celerity
StationTime(min)(deg)(m/s)(Hz)(m/s)(m/s)
I48TN 17:06:41 21:13.1 81.5–93.0 336.0–392.0 0.2–6.4 285.3–336.4 297 
I26DE 17:12:21 15:25.8 119.2–134.4 328.0–401.0 0.1–2.0 292.7–329.1 297 
I17CI 19:44:32 13:33.1 44.0–50.0 302.0–364.0 0.1–1.6 293.6–308.1 314 

The modeling of the infrasound propagation from the source to receiver in this paper uses normal modes and the PE with the range-independent and range-dependent models. The ground to space vertical wind profile is the main input for the infrasound propagation simulation. Furthermore, the temperature, density, molecular mass, and pressure as a function of the altitude (z) are also needed for the modeling. The Modern-Era Retrospective Analysis for Research and Applications version 2 (MERRA-2) reanalysis data are used in both the semiempirical and parabolic transmission loss estimations. The reconstruction of the wind profile depends on the spline interpolation and estimation of the spherical harmonics (SH), and the decomposition between the different wind data sets by the AVOG2S model (Schwaiger et al., 2019) were calculated. The wind data for higher altitudes (more than 76 km) is extracted from the Horizontal Wind Model 14 (HWM-14; Drob et al., 2015). The profile of the temperature, density, molecular mass, and pressure are estimated from the empirical model of the MSISE (Picone et al., 2002). The ground to space (G2S; Drob et al., 2003) vertical wind profiles [Fig. 6(a)] above the event location show some jets to the north-east in the troposphere up to 20 km. In contrast, the wind jets in the stratosphere are strongly to the west, which suggests effective infrasound ducting in this direction. Figure 6(b) compares the dominant wind speed and direction obtained from the reanalysis of the European Centre for Medium-Range Weather Forecasts (Hersbach et al., 2018) at an altitude of 48 km, which reflects the strong westward wind jets with a speed near 41 m/s over the B-p.

FIG. 6.

(Color online) The (a) extracted G2S vertical wind profile above the Beirut explosion site and (b) dominant wind speed over the region at an altitude of 48 km, derived from European Centre for Medium-Range Weather Forecasts (ECMWF) at 15:00:00 of 4 August 2020, are depicted.

FIG. 6.

(Color online) The (a) extracted G2S vertical wind profile above the Beirut explosion site and (b) dominant wind speed over the region at an altitude of 48 km, derived from European Centre for Medium-Range Weather Forecasts (ECMWF) at 15:00:00 of 4 August 2020, are depicted.

Close modal

The thermodynamic speed cth was estimated from the equation (z)=γRT (Garcés et al., 1998a), where γ (ratio of specific heats) = 1.4 and R (the specific gas constant for dry air) = 286.9 J kg−1 K−1. Moreover, the effective sound speed ceff(z) can be obtained from the relation ceff(z)=cth(z)+|wuv(z)cos(ϕϕwuv(z))|, where wuv is the horizontal wind, ϕ is the azimuth direction of propagation, and ϕwuv is the wind direction (Smets et al., 2016). Figure 7 shows the estimation of ceff(z) in all directions from the source for every 45 deg.

FIG. 7.

(Color online) The calculated effective sound speed (red lines) and thermodynamic sound speed (blue lines) are shown. The brown vertical lines refer to the adiabatic sound speed at zero altitude (340 m/s).

FIG. 7.

(Color online) The calculated effective sound speed (red lines) and thermodynamic sound speed (blue lines) are shown. The brown vertical lines refer to the adiabatic sound speed at zero altitude (340 m/s).

Close modal

The normal mode expansion method with an approximation of the effective sound speed (Waxler et al., 2017a) was applied using the frequency of 0.5 Hz along all azimuths from 0 to 360 deg to a maximum distance of 6000 km to estimate the power transmission loss of the infrasound waves in all directions as shown in Fig. 8. This method is effective for the low angle propagation paths and low wind speeds; moreover, the accurate results for the tropospheric and stratospheric ducting can be given from the perturbative treatment of the atmospheric attenuation. However, the accuracy of the returned signal from the thermosphere is not as robust.

FIG. 8.

(Color online) The 2D-transmission loss for 0.5 Hz estimate by the normal modes (range-independent) method in all directions is shown. Maps Data: Google, Data SIO, NOAA, U.S. Navy, NGA, GEBCO; Image NOAA; Image Landsat; Data USGS.

FIG. 8.

(Color online) The 2D-transmission loss for 0.5 Hz estimate by the normal modes (range-independent) method in all directions is shown. Maps Data: Google, Data SIO, NOAA, U.S. Navy, NGA, GEBCO; Image NOAA; Image Landsat; Data USGS.

Close modal

To solve the Helmholtz equation, we let ρ0 be the mean density of the atmosphere, α is the atmospheric attenuation coefficient, ω is the angular frequency, and p is the deviation from the mean pressure at the horizontal position xH and altitude z,

[H2+ρ0z(1ρ0z)+(ωceff(z)+iα(ω,z))2]p(r,z,ω)=0.
(5)

p(r,z,ω) can be estimated from Eq. (5), depending on the variation in the effective sound speed and the azimuth direction of the estimation. Therefore, the transmission loss can be calculated as

I(r,z,ω)=18πrρ0(z)ρ0(zs)je2αjrkj(ψj(zs)ψj(z))2.
(6)

A stratospheric duct exists to the west from the event site with a transmission loss of 0.5 Hz in the range from 75 to 101 dB, suggesting that the infrasound waves might be observed at the HLW, I26DE, I48TN, I42PT, I17CI, and I11CV infrasound stations. However, the transmission loss effectively increases in the red and orange zones in Fig. 7, indicating that the infrasound waves in these zones might not be recorded. The range-dependent PE method is also used to estimate the two-dimensional (2D)-transmission loss for 0.5 Hz in the directions of the I48TN, I26DE, and I17CI IMS (International Monitoring System) infrasound stations as shown in Fig. 9. Transmission losses of <130 dB are predicted for I48TN and I26DE, but in the case of I17CI, the transmission loss of >150 dB. The thermospheric phase is expected to disappear at about 4800 km as a result of the absorption and long path effect.

FIG. 9.

(Color online) The 2D-transmission loss simulations by the parabolic range-dependent propagation method in the directions of (a) I48TN, (b) I26DE, and (c) I17CI are shown.

FIG. 9.

(Color online) The 2D-transmission loss simulations by the parabolic range-dependent propagation method in the directions of (a) I48TN, (b) I26DE, and (c) I17CI are shown.

Close modal

According to Zhang et al. (2019), at 14:48 China local time (UTC+8), a first fire flame was observed and followed by a small explosion. A few seconds later, a loud explosion was heard for several seconds and the flame heights reached 30 m. In addition, the buildings within 500 m were totally destroyed. Furthermore, damage of numerous houses was reported within 2 km of the explosion site. According to Ramzy and Hernández (2019), within 5 km of the explosion, many of the windows of the houses were broken.

The Kochi University of Technology (KUT) sensor network consists of 30 infrasound sensors, which were established as part of the tsunami early warning system in Shikoku Island (Hamama and Yamamoto, 2021; Saito et al., 2021). The explosion was well-recorded on most of these sensors especially on Shikoku Island (K11, K12, K13, K14, K16, K17, and K41) as shown in Fig. 10.

FIG. 10.

(Color online) The infrasound signals observed at K11, K12, K13, K14, K16, K17, and K41 on the KUT sensor network for the Xi chemical explosion are shown.

FIG. 10.

(Color online) The infrasound signals observed at K11, K12, K13, K14, K16, K17, and K41 on the KUT sensor network for the Xi chemical explosion are shown.

Close modal

The computed celerties for the infrasound arrivals of the KUT sensors are about 338 m/s, which confirms ducting in the tropospheric layer as shown in Fig. 11.

FIG. 11.

(Color online) The travel-time curve between the Xi chemical explosion and KUT sensors network, where the blue line refers to a celerity of 338 m/s, is shown.

FIG. 11.

(Color online) The travel-time curve between the Xi chemical explosion and KUT sensors network, where the blue line refers to a celerity of 338 m/s, is shown.

Close modal

The G2S profile over the explosion site was extracted from the reanalysis of MERRA-2 for the case of Xi. In contrast to the B-p explosion, strong tropospheric jets were clearly observed to the east, where the zonal wind speed reached a maximum of 63 m/s at an altitude of 13.4 km. These strong jets can strengthen the infrasound propagation to the direction of the eastern side. Normal mode modeling was applied to estimate the transmission losses at a 0.5 Hz frequency along all azimuths between 0 and 360 deg to a distance of 2000 km as shown in Fig. 12, which confirmed the propagation to the east of the explosion as the transmission loss-estimates were less than 100 dB. Furthermore, the PE modeling was investigated along an azimuth of 94 deg, consistent with the KUT sensor network (Fig. 13), confirming the arrival of the tropospheric phases to the sites in Japan.

FIG. 12.

(Color online) The 2D-transmission loss for 0.5 Hz from the Xi chemical explosion site is shown. Maps Data: Google, Data SIO, NOAA, U.S. Navy, NGA, GEBCO; Image NOAA; Image Landsat; Data USGS.

FIG. 12.

(Color online) The 2D-transmission loss for 0.5 Hz from the Xi chemical explosion site is shown. Maps Data: Google, Data SIO, NOAA, U.S. Navy, NGA, GEBCO; Image NOAA; Image Landsat; Data USGS.

Close modal
FIG. 13.

(Color online) The 2D-transmission loss-estimates using a range-dependent parabolic-equation propagation model in the directions of the KUT sensor network are shown. Moreover the G2S profile shows strong tropospheric wind to the east direction.

FIG. 13.

(Color online) The 2D-transmission loss-estimates using a range-dependent parabolic-equation propagation model in the directions of the KUT sensor network are shown. Moreover the G2S profile shows strong tropospheric wind to the east direction.

Close modal

The yield estimation from the infrasound has a long history. Pierce and Posey (1971) investigated the use of multiple waveform measures, including the zero-to-peak pressure (P) and time of the first complete cycle of the waveform (T) with the angular distance Δ with the yield (W) estimated using the relation: W=0.494P(sinΔ)1/2Hs(T)3/2, where Hs is the atmospheric scale height. The relation was applied to the data from the atmospheric nuclear tests and found to provide good estimates of the yield. In addition, the source duration (Ts) can be estimated from the relation Ts=0.33W1/3 s Pierce and Posey (1971), dependent on using the first cycle of the Lamb wave, which can be defined as a surface guided wave that can accompany tropospheric arrivals from the explosive sources (Garcés et al., 1998; ReVelle and Kulichkov, 1998). Furthermore, modification of their relationship was documented by Blandford and Clauter (1995) and is known as the AFTAC relationship. They similarly assumed that the period (T) is independent of the distance and proportional to the cubic root of the yield. Moreover, the absolute pressure is proportional to the square root of the distance within 20 deg (Stevens et al., 2002).

According to Green and Bowers (2010), both the AFTAC relation and Whitaker (1995), which is known as the LANL relationship, can be applied to the chemical and nuclear explosions with yields less than 2 kt. Furthermore, several experiments conducted had evaluated the yield estimates from the seismo-acoustic signals for the near-surface chemical explosions, and the charge depends on the height of the burst or depth of the burial for the explosion through the inversions of the scaled seismic P wave displacement and the air overpressure (Ford et al., 2014).

In this study, the Whitaker (1995) relationship, Pc=2.35.103(R/W1/2)1.36, was used for the infrasound yield estimation. The advantage of the LANL relationship is that several experiments were performed with a variety of noise models, different infrasound propagation, as well as different noise reduction systems with the attenuation considered above 10 km. In addition, the wind speed has a significant effect on the sound speed at these heights (>10 km), hence, it can affect the propagation speed and surface bounce locations. Thus, the atmospheric wind is a part of the yield estimation (Stevens et al., 2002; Whitaker, 1995). Pc is the corrected zero-to-peak pressure and R is the distance in kilometers. Here, a correction based on the stratospheric wind speed was also taken into account for the observed zero-to-peak pressure P using Pc=100.019vP, where Pc is the corrected zero-to-peak pressure and v is the stratospheric wind speed to the direction of the wave motion. The stratospheric wind speed to the west at the time of the B-p explosion was 48 m/s, and for the Xi explosion, the stratospheric wind speed in the direction of the KUT sensor network was 32 m/s. In addition, according to the relation of the International Data Center (IDC) by Brown (1999), the infrasound magnitude was subtracted from the estimated yield to the following relations (Stevens et al., 2006):

Corr.Mag.=log(P)+1.36log(R)0.019v,

where W is the charge in kt, P is the max zero-to-peak pressure, R is the distance in kilometers, and v is the stratospheric wind speed at an altitude of 50 km. The previous relation is originally derived from the LANL relation,

Mag.=0.68log(W)+3.37,

where the magnitude term can represent the source size with independence of the distance at which the signal is measured. The mean yield estimate of the B-p explosion using the data from the IMS station is 271 (± 101.12) tons (Table II), whereas the yield estimate for the Xi explosion is 201.6 (± 83.51) tons (Table III).

TABLE II.

The estimated yield from the recorded infarsonic waves of the B-p explosion.

DistanceObserved maximun pressureEstimated yield
Station(km)(Pa)(ton)IDC magnitudeIDC magnitude corrected
I48TN 2390 0.188 230.6 2.372 2.9379 
I26DE 2450 0.143 172.8 2.1762 2.852 
I17CI 5100 0.095 410.3 2.76 3.108 
DistanceObserved maximun pressureEstimated yield
Station(km)(Pa)(ton)IDC magnitudeIDC magnitude corrected
I48TN 2390 0.188 230.6 2.372 2.9379 
I26DE 2450 0.143 172.8 2.1762 2.852 
I17CI 5100 0.095 410.3 2.76 3.108 
TABLE III.

The estimated yield from the recorded infarsonic waves of the Xi chemical explosion.

DistanceObserved maximum pressureEstimated yield
Station(km)(Pa)(ton)MagnitudeMagnitude corrected
K12 1270 0.100 107.0 1.50 2.55 
K41 1118 0.3100 309.0 2.57 3.03 
K13 1278 0.158 152.0 2.09 2.81 
K11 1256 0.341 295.0 2.83 3.13 
K14 1308 0.148 145.0 2.05 2.80 
DistanceObserved maximum pressureEstimated yield
Station(km)(Pa)(ton)MagnitudeMagnitude corrected
K12 1270 0.100 107.0 1.50 2.55 
K41 1118 0.3100 309.0 2.57 3.03 
K13 1278 0.158 152.0 2.09 2.81 
K11 1256 0.341 295.0 2.83 3.13 
K14 1308 0.148 145.0 2.05 2.80 

The satellite based remote sensing technologies using Interferometric Synthetic Aperture Radar (InSAR) coherence to calculate the surface deformation were applied to the explosions in this study. The InSAR methodology depends on the measurement of the phase difference between two interferograms acquired at different times. In recent decades, this technique was proven to be useful in estimating the surface land subsidence and uplifts that can reach a few millimeters of displacement (Massonnet et al., 1993; Wempen, 2020). The InSAR coherence has been used to map the environmental damages from landslides, volcanic eruptions, and burnt scars and their severity (Gaebler et al., 2019; Spaans and Hooper, 2016). The European Space Agency's (ESA) Sentinel-2 optical data with 10 m spatial resolution have been used to rapidly map the effects of the fire, urbanization, the normalized difference vegetation index (NDVI), and land subsidence (e.g., Martinis et al., 2017; Tzouvaras et al., 2020; Xu et al., 2018).

The synthetic aperture radar (SAR) C-band and the chromatic bands with high spatial and temporal resolution, which were acquired from Sentinel-1/2 (S-1 and S-2) satellite sensors (Attema et al., 2010) are used to determine the surface deformation accompanying the explosions. Sentinel-1 radar data S-1 and Sentinel-2 Multispectral Imager (MSI) S-2 provided 5.6 cm resolution and 5.4 GHz frequency for the radar scenes. Thirteen optical bands in the visible, near-infrared (NIR), and shortwave infrared (SWIR) spectrum were used. Outputs from the two sensors were compared before and after the event to map the areas of deformation, where the explosion desolated most. The InSAR coherence maps for the B-p used the data collected between 18 July and 11 August 2020 from the C-band Sentinel-1 pair images in the Single Look Complex (SLC) format acquired in the Interferometric Wide Swath (IW) TOPSAR mode from the descending satellite tracks. The interferograms were selected with small perpendiculars to the baseline orbit of less than 150 m and threshold coherence of 0.2.

The Sentinel-1 interferograms were registered precisely using the enhanced spectral diversity (ESD) algorithm to mitigate the phase discontinuities resulting from the ionospheric propagation effects.The interferograms were generated after interferometric phase noise reduction, including flat earth and topographic phase removal, followed by multi-looking processing techniques to enhance the amplitude quality. A window size of 5 × 5 m was used to improve the coherency by spatial averaging. The range Doppler terrain correction (RDTC) was applied to the stacked interferograms. The calculated interferometric coherence between the two SLC look images is defined as the complex amplitude of the correlation coefficient between the master and slave scenes. The interferometric coherence value from 0 and 1, as shown in Eq. (7), is

γ=|γ|ejϕ=|S1S2|S1S1S2S2,
(7)

where S1 and S2 are the complex SLC looks of the master and slave scenes, S* is its conjugate. To determine the degree of surface change due to the damage, the normalized interferometric coherence difference was used as represented in Eq. (8),

γDiff=γpreγafterγpre+γafter,
(8)

where γpre, γafter, and γDiff are the coherence parameters of the SAR pair images before and after the explosion, and the normalized coherence difference (NCD) between γpre and γafter, respectively (Molan and Lu, 2020). The coherence change detection (CCD) method derived from the NCD proved to be more effective than the coherence difference (CD) in estimating the areas affected by the damage through minimizing the false damage pixels (see Tzouvaras et al., 2020). According to the normalized difference of the coherency change value, the explosion damage was classified into three categories, the highest level of damage around the blast site area ranged between 1 and 0.6 and is represented by heavy damage, and the medium damage with the NCD values between 0.6 and 0.2 were annotated with the red to orange colors in Fig. 14(b) and the lowest NCD values (less than 0.2) refer to slight and none-damage. In addition, these results were confirmed by the correlation with ground truth data observed in this area, updated information reported on Google Earth (computer program, Google, Mountain View, CA), and available news sources. Finally, the building location-maintained values, which were lower than 0.2 to -1, corresponded to slight damage, in which the urban buildings in these areas reported glass window damage and minor distortion as shown in Fig. 14 (Tamkuan and Nagai, 2017; Watanabe et al., 2016). Because the study area is surrounded by the Mediterranean Sea, pixels from the water are removed. However, the Normalized Difference Water Index (NDWI) masked out the effect of the water pixels (Efthimiou et al., 2020; Fig. 14). Three deformation areas had been reported: the confirmed total infrastructure damage area, extensive infrastructure damage area, and reported minor damage area. Land deformation reached 20 km from the blast site. The surface destruction estimated from Sentinel-1/2 confirmed the total infrastructure damage reported area, which extended 3 km, concentrated to the southwest direction of the blast site.

FIG. 14.

(Color online) Examples of the damaged areas of the (a) Xi and (b) B-p cases before and after the explosion. The ground truth confirmed the estimated damage level, and the accuracy of the ground truth was observed from Google Earth (Google, Mountain View, CA) and the reported news. Maps Source: Esri, Maxar, GeoEye, Earthstar Geographics, CNES/Airbus DS, USDA, USGS, AeroGRID, IGN, and the GIS User Community.

FIG. 14.

(Color online) Examples of the damaged areas of the (a) Xi and (b) B-p cases before and after the explosion. The ground truth confirmed the estimated damage level, and the accuracy of the ground truth was observed from Google Earth (Google, Mountain View, CA) and the reported news. Maps Source: Esri, Maxar, GeoEye, Earthstar Geographics, CNES/Airbus DS, USDA, USGS, AeroGRID, IGN, and the GIS User Community.

Close modal

We also applied the previous processing steps to the Sentinel-1 images to determine the impacts of the Xi chemical explosion. The interferometric coherence map was constructed from the interferometric pairs dated 7 March and 31 March 2019. According to the location of the blast site shown in Fig. 14(b), the site is surrounded by a high vegetation density area that distorts the Sentinel-1 back scattering wave; therefore, the NDVI obtained from the Sentinel-2 optical satellite (level-2A) images, acquired on 15 March 2019 (Braun et al., 2019), followed by the NDVI mask were used as the calculation to remove the coherence loss effect (NDVI 0.2) from the output map (Fig. 14).

To assess the comparison between the NCD observation and its correspondence to the explosion yield for the B-p and Xi chemical explosion damages, the empirical boom relation was proposed. Figures 19 and 20 (see the  Appendix) illustrate the differences between the damages produced by the two explosions, in which the Xi explosion showed high NCD values corresponding to the destructive damage in the closest area 600 m from the explosion. In contrast, the B-p explosion warehouse site was close to the sea and based on the location conditions, the impact of the B-p explosion was less, showing a high NCD value to a distance of 300 m away from the explosion.

The peak overpressure based on the explosion yield can be estimated using the Boom with the damage term estimated through the relation between the NCD and peak overpressure (P) in kilopascals with the distance (r) in meters as shown in the Eq. (9), where the damage proxy (D) can be estimated from Eq. (10) (Pilger et al., 2021),

P=3.45978103W0.444A0.556r1.333,
(9)
D=log10P(r)100log10max(P),
(10)

where A is the atmospheric pressure and W is the charge in kilotons. We applied Eqs. (9) and (10) to both explosions, considering that the atmospheric temperature effect was not taken into account at each point of the peak overpressure estimation. We assumed the maximum peak overpressure in the damage relation to be 80 kPa in the B-p explosion and 88 kPa in Xi explosion, according to the dominant maximum overpressure, which occurred in the damaged area with the highest NCD percentage (NCD%). Figure 15 displays the relation between the damage and calculated overpressure for both explosions with the yields ranging from 0.1 to 2 kt in the B-p explosion. The relation is well-correlated with the yield range of 0.1–0.6 kt in the area 600 m from the explosion although the logarithmic trend of the NCD% follows the higher charge estimates (0.6–1 kt). However, in the Xi case, the damage proxy located in the area of 300 m does not follow the NCD% trend, although after 300 m to 2 km, the NCD% logarithmic trend is consistent with a charge of 0.1–1 kt.

FIG. 15.

(Color online) The relationships between the NCD% and damage percentage in the case of the B-p explosion (a) and Xi explosion (b) using a charge-range from 0.1 to 2 kt are shown.

FIG. 15.

(Color online) The relationships between the NCD% and damage percentage in the case of the B-p explosion (a) and Xi explosion (b) using a charge-range from 0.1 to 2 kt are shown.

Close modal

In the case of the Beirut explosion, the open access seismic stations and Egyptian National Seismic Network have recorded the acoustic waves to distances up to 316 km, where the seismic and hydroacoustic signals are observed to distances of up to 727.1 km. The estimated magnitude from the 20 stations is 3.55 (± 0.15 ML). In addition, the infrasound propagation confirms a stratospheric duct toward I48TN, I26DE, and I17CI. Moreover, the estimated explosion equivalent TNT yield from the seismic observations is 202.2 (± 127.55 tons). Furthermore, the infrasound yield estimation equivalent to TNT is 271 (± 100 tons). The analysis of the remote sensing data shows three deformation areas: the confirmed total infrastructure damage, extensive infrastructure damage, and reported minor damage. The surface destruction estimated from Sentinel-1/2 confirms the total of the reported infrastructure damage area, which extends 2 km. In addition, the estimated yield range was evaluated from the NCD% pixels and modeled overpressure values to 2 km from the explosion site. The estimated damage from the overpressure values and the logarithmic trend of the NCD% are well-fit with some limitation to the yield range of 0.1–0.6 kt out to the distance of 600 m, and after that, the NCD% logarithmic trend follows the higher yield estimates (0.6–1 kt).

Converting the yield into an equivalent weight, ammonium nitrate is more complicated. The ammonium nitrate exists in several product categories of different densities and concentrations and is also being filled in different package sizes. For the typical low-density ammonium nitrate, in stockpiles with an explosion due to fire, the overall TNT equivalency value according to GHD (2012) is 10%, 15%, and 20% in bulk in the free stockpiles, bays, and bags, respectively.

The ammonium nitrate equivalent yield average (from the seismic and infrasound yield estimations) is between 2020 and 2700 tons. The declared amount stored in the port was around 2700 tons of ammonium nitrate. The difference between the estimated and declared ammonium nitrate masses could have been burned during the fire before the explosion. Otherwise, the ammonium nitrate could have been hydrated during the long period of storage in poor conditions thereby reducing the W explosive yield.

In the Xi chemical explosion case, open access seismic data sets were not available. The strong eastward tropospheric wind proved well-recorded infrasound waves at the KUT sensor network in Japan. The infrasound yield estimate equivalent to TNT is 201 (± 83.5 tons). The estimated yield range was evaluated using the InSAR normalized interferometric coherence pixels and modeled overpressure values to a 2 km range from the explosion site. The estimated damage from the overpressure values and logarithmic trend of the NCD% are in the range of yield of 0.1–1 kt after a distance of 300 m. This relation needs improvements to assess the building distributions and conditions of the damage at each site. In addition, further limitations can be overcome by using the multi-temporal coherence analysis with a longer temporal baseline (days). This study, using multiple techniques, quantifies the effects of the near-surface explosions and safety area from the chemical factory or blasts needed to avoid the loss of life and destruction of buildings.

The wind datasets were extracted from Modern-Era Retrospective analysis for Research and Applications, version 2 (MERRA-2), which can be accessed online).1 Some of the figures were created with Generic Mapping Tools (GMT) and Quantum Geographic Information System (QGIS) 3.6 mapping programs with the base maps from Google Maps (Google, Mountain View, CA). The seismic data were obtained from the Egyptian National Seismic Network (ENSN) and the GEOFON Data Centre of the GeoForschungsZentrum (GFZ) German Research Centre for Geosciences. We acknowledge the forthcoming Sentinel-1 satellite designed by the European Space Agency (ESA). In addition, the infrasound data were extracted from the Kochi University of Technology (KUT) Sensor Network, which is accessible online.2 Furthermore, the International Monitoring System (IMS) infrasound data sets of Comprehensive Nuclear Test Ban Treaty (CTBTO) were extracted through the secure web portal of CTBTO.3 The authors would like to thank the anonymous reviewers and Dr. D. Keith Wilson for their valuable comments, which helped significantly in improving the content of the manuscript. The researchers I.H. and N.I.M. are funded by the Egypt-Japan Education Partnership “EJEP-3” scholarship from the Ministry of Higher Education of the Arab Republic of Egypt. This work was partly supported by the SECOM Science and Technology Foundation and JSPS KAKENHI Grant No. JP17H02062.

B-p

Beirut-Port

Xi

Xiangshui

max

Maximum

avg

Average

Tables IV and V depict the of seismic and acoustics arrivals and the local magnitude estimations. Figure 16 show the waveform and spectrogram for signals from some seismic stations. Moreover, Figs. 17 and 18 refer to the estimated location from seismic and infrasound stations. In addition, Figs. 19 and 20 show examples of ground truth data images of B-p and Xi explosions.

TABLE IV.

The arrivals of the recorded seismic phases at stations from the Egyptian National Seismic Network (HL) and other open access networks (IS,IM,CQ,GE,andKO).

DistanceP-phase arrivalsT-phase arrivalsShock wave arrivals
NumberStationNetwork(km)(UTC)(UTC)(UTC)
GEM IS 77.70 15:08:33 No phase No phase 
CY606 IM 103.23 15:08:37 15:09:26 No phase 
CY602 IM 104.34 15:08:38 15:09:27 No phase 
CY604 IM 104.34 15:08:39 15:09:27 No phase 
CY603 IM 105.45 15:08:39 15:09:27 No phase 
CY601 IM 105.45 15:08:39 15:09:27 No phase 
CY401 IM 148.74 15:08:43 15:09:55 No phase 
MMLIS 163.17 15:08:47 No phase No phase 
PARA CQ 182.24 15:08:47 15:10:16 No phase 
10 SALP GE 205.35 15:08:53 No phase No phase 
11 MVOU CQ 212.97 15:08:53 15:10:35 15:20:18 
12 UJAP GE 216.45 15:08:53 No phase No phase 
13 EREN KO 218.67 15:08:53 15:10:38 15:20:42 
14 ASGA CQ 229.77 15:08:54 15:10:41 15:21:09 
15 CSS GE 232.85 15:08:54 15:10:37 15:21:16 
16 CY303 IM 237.54 15:08:54 15:10:51 No phase 
17 LFK KO 237.54 15:08:54 15:10:41 15:21:31 
18 ATHA CQ 237.54 15:08:55 15:11:30 15:21:33 
19 CY302 IM 238.65 15:08:55 15:10:52 No phase 
20 DSI IS 256.41 15:08:59 No phase No phase 
21 AMAS IS 268.62 15:09:01 No phase No phase 
22 THAT KO 280.83 15:09:02 No phase No phase 
23 YTIR IS 284.16 15:09:03 No phase No phase 
24 NATA CQ 287.49 15:09:03 No phase 15:23:52 
25 MSBI GE 287.49 15:09:03 No phase No phase 
26 GHAJ GE 288.76 15:09:04 No phase No phase 
27 ALEF CQ 301.78 15:09:05 15:11:55 15:24:32 
28 AKMS CQ 316.35 15:09:06 No phase 15:25:16 
29 CY201 IM 334.11 15:09:07 15:11:57 No phase 
30 KZIT IS 348.54 15:09:10 No phase No phase 
31 CY101 IM 359.64 15:09:10 15:12:14 No phase 
32 PRNI IS 397.38 15:09:16 No phase No phase 
33 KRMI IS 426.24 15:09:20 No phase No phase 
34 HRFI IS 430.68 15:09:21 No phase No phase 
35 EIL GE 472.86 15:09:25 No phase No phase 
36 URFA KO 492.84 15:09:27 No phase No phase 
37 SUZ HL 517.86 15:09:31 No phase No phase 
38 BST HL 526.22 15:09:31 No phase No phase 
39 NUB HL 552.53 15:09:36 No phase No phase 
40 KOT HL 562.51 15:09:36 No phase No phase 
41 DHB HL 582.22 15:09:39 No phase No phase 
42 KAT HL 615.39 15:09:42 No phase No phase 
43 NDB3 HL 727.13 15:09:55 15:17:10 No phase 
DistanceP-phase arrivalsT-phase arrivalsShock wave arrivals
NumberStationNetwork(km)(UTC)(UTC)(UTC)
GEM IS 77.70 15:08:33 No phase No phase 
CY606 IM 103.23 15:08:37 15:09:26 No phase 
CY602 IM 104.34 15:08:38 15:09:27 No phase 
CY604 IM 104.34 15:08:39 15:09:27 No phase 
CY603 IM 105.45 15:08:39 15:09:27 No phase 
CY601 IM 105.45 15:08:39 15:09:27 No phase 
CY401 IM 148.74 15:08:43 15:09:55 No phase 
MMLIS 163.17 15:08:47 No phase No phase 
PARA CQ 182.24 15:08:47 15:10:16 No phase 
10 SALP GE 205.35 15:08:53 No phase No phase 
11 MVOU CQ 212.97 15:08:53 15:10:35 15:20:18 
12 UJAP GE 216.45 15:08:53 No phase No phase 
13 EREN KO 218.67 15:08:53 15:10:38 15:20:42 
14 ASGA CQ 229.77 15:08:54 15:10:41 15:21:09 
15 CSS GE 232.85 15:08:54 15:10:37 15:21:16 
16 CY303 IM 237.54 15:08:54 15:10:51 No phase 
17 LFK KO 237.54 15:08:54 15:10:41 15:21:31 
18 ATHA CQ 237.54 15:08:55 15:11:30 15:21:33 
19 CY302 IM 238.65 15:08:55 15:10:52 No phase 
20 DSI IS 256.41 15:08:59 No phase No phase 
21 AMAS IS 268.62 15:09:01 No phase No phase 
22 THAT KO 280.83 15:09:02 No phase No phase 
23 YTIR IS 284.16 15:09:03 No phase No phase 
24 NATA CQ 287.49 15:09:03 No phase 15:23:52 
25 MSBI GE 287.49 15:09:03 No phase No phase 
26 GHAJ GE 288.76 15:09:04 No phase No phase 
27 ALEF CQ 301.78 15:09:05 15:11:55 15:24:32 
28 AKMS CQ 316.35 15:09:06 No phase 15:25:16 
29 CY201 IM 334.11 15:09:07 15:11:57 No phase 
30 KZIT IS 348.54 15:09:10 No phase No phase 
31 CY101 IM 359.64 15:09:10 15:12:14 No phase 
32 PRNI IS 397.38 15:09:16 No phase No phase 
33 KRMI IS 426.24 15:09:20 No phase No phase 
34 HRFI IS 430.68 15:09:21 No phase No phase 
35 EIL GE 472.86 15:09:25 No phase No phase 
36 URFA KO 492.84 15:09:27 No phase No phase 
37 SUZ HL 517.86 15:09:31 No phase No phase 
38 BST HL 526.22 15:09:31 No phase No phase 
39 NUB HL 552.53 15:09:36 No phase No phase 
40 KOT HL 562.51 15:09:36 No phase No phase 
41 DHB HL 582.22 15:09:39 No phase No phase 
42 KAT HL 615.39 15:09:42 No phase No phase 
43 NDB3 HL 727.13 15:09:55 15:17:10 No phase 
TABLE V.

The estimation of the local magnitude from the three components of the recorded seismic waveforms.

A(N)A(E)DistantceMLMLML
Station(mm)(mm)(km)(N component)(E component)(mean)log(A0)
CY606 2.196 2.774 103 3.361 3.463 3.412 3.020 
CY602 2.785 1.998 104 3.471 3.327 3.399 3.026 
CY604 1.887 1.778 104 3.302 3.276 3.289 3.026 
CY603 3.124 2.272 106 3.534 3.396 3.465 3.039 
CY601 3.139 2.149 106 3.536 3.372 3.454 3.039 
MMLI 1.366 1.119 163 3.490 3.404 3.447 3.355 
SALP 1.222 — 205 3.631 — 3.631 3.544 
UJAP 1.210 1.094 216 3.673 3.629 3.651 3.590 
CY303 0.566 0.499 237 3.428 3.373 3.400 3.675 
CY302 0.806 0.608 239 3.589 3.467 3.528 3.683 
AMAZ 0.677 0.689 269 3.627 3.635 3.631 3.796 
YTIR 0.682 0.407 284 3.685 3.460 3.573 3.851 
GHAJ 1.679 1.002 289 4.094 3.870 3.982 3.869 
CY201 0.458 0.549 334 3.684 3.763 3.724 4.024 
KZIT 0.258 0.364 348 3.481 3.630 3.556 4.070 
KRMI 0.191 0.203 426 3.596 3.621 3.609 4.315 
HRFI 0.168 0.202 431 3.556 3.634 3.595 4.330 
EIL 0.074 0.079 473 3.324 3.352 3.338 4.454 
URFA 0.140 0.119 493 3.658 3.586 3.622 4.512 
BST 0.142 0.114 527 3.761 3.666 3.714 4.608 
A(N)A(E)DistantceMLMLML
Station(mm)(mm)(km)(N component)(E component)(mean)log(A0)
CY606 2.196 2.774 103 3.361 3.463 3.412 3.020 
CY602 2.785 1.998 104 3.471 3.327 3.399 3.026 
CY604 1.887 1.778 104 3.302 3.276 3.289 3.026 
CY603 3.124 2.272 106 3.534 3.396 3.465 3.039 
CY601 3.139 2.149 106 3.536 3.372 3.454 3.039 
MMLI 1.366 1.119 163 3.490 3.404 3.447 3.355 
SALP 1.222 — 205 3.631 — 3.631 3.544 
UJAP 1.210 1.094 216 3.673 3.629 3.651 3.590 
CY303 0.566 0.499 237 3.428 3.373 3.400 3.675 
CY302 0.806 0.608 239 3.589 3.467 3.528 3.683 
AMAZ 0.677 0.689 269 3.627 3.635 3.631 3.796 
YTIR 0.682 0.407 284 3.685 3.460 3.573 3.851 
GHAJ 1.679 1.002 289 4.094 3.870 3.982 3.869 
CY201 0.458 0.549 334 3.684 3.763 3.724 4.024 
KZIT 0.258 0.364 348 3.481 3.630 3.556 4.070 
KRMI 0.191 0.203 426 3.596 3.621 3.609 4.315 
HRFI 0.168 0.202 431 3.556 3.634 3.595 4.330 
EIL 0.074 0.079 473 3.324 3.352 3.338 4.454 
URFA 0.140 0.119 493 3.658 3.586 3.622 4.512 
BST 0.142 0.114 527 3.761 3.666 3.714 4.608 
FIG. 16.

(Color online) The spectrogram for the recorded signal to calculate the frequency content is shown.

FIG. 16.

(Color online) The spectrogram for the recorded signal to calculate the frequency content is shown.

Close modal
FIG. 17.

(Color online) The detected seismic stations of the B-p explosion are shown.

FIG. 17.

(Color online) The detected seismic stations of the B-p explosion are shown.

Close modal
FIG. 18.

(Color online) The cross-bearing is shown by plotting the azimuths between the sites. Maps Data: Google, Data SIO, NOAA, U.S. Navy, NGA, GEBCO; Image Landsat / Copernicus; Image IBCAO.

FIG. 18.

(Color online) The cross-bearing is shown by plotting the azimuths between the sites. Maps Data: Google, Data SIO, NOAA, U.S. Navy, NGA, GEBCO; Image Landsat / Copernicus; Image IBCAO.

Close modal
FIG. 19.

(Color online) Examples of the damaged areas before (left column) and after (right column) the Beirut explosion. The estimated damage level is confirmed by the ground truth images. The A, B and C panels (corresponding to each damage category on the map shown in Fig. 14), extracted from Google Earth (Google, Mountain View, CA), emphasized with the degree of damage, correspond to the values observed on the interferometric coherence maps; the destructive damage ranges from 1 to 0.6, and the heavy to slight damage of interferometric coherence range from 0.6 to 0.2 located in the area surrounding the source of the explosion and are annotated with labels (B and C). Maps Data: (A–B) Google, image © 2021 CNES /Airbus; (C) Google, image © 2021 Maxar Technologies.

FIG. 19.

(Color online) Examples of the damaged areas before (left column) and after (right column) the Beirut explosion. The estimated damage level is confirmed by the ground truth images. The A, B and C panels (corresponding to each damage category on the map shown in Fig. 14), extracted from Google Earth (Google, Mountain View, CA), emphasized with the degree of damage, correspond to the values observed on the interferometric coherence maps; the destructive damage ranges from 1 to 0.6, and the heavy to slight damage of interferometric coherence range from 0.6 to 0.2 located in the area surrounding the source of the explosion and are annotated with labels (B and C). Maps Data: (A–B) Google, image © 2021 CNES /Airbus; (C) Google, image © 2021 Maxar Technologies.

Close modal
FIG. 20.

(Color online) The Xi ground truth extracted from Google Earth (Google, Mountain View, CA) before and after the explosion (corresponding to each symbol on the map shown in Fig. 14). The destructive damage annotated with (A) has the highest coherence value range from 0.6 to 1, whereas (B) and (D) explain the rate of damage according to the distance to the explosion. The slight damage is annotated with (D) showing the lowest interferometric coherence value 0.2. Maps Data: Google, image © 2021 Maxar Technologies.

FIG. 20.

(Color online) The Xi ground truth extracted from Google Earth (Google, Mountain View, CA) before and after the explosion (corresponding to each symbol on the map shown in Fig. 14). The destructive damage annotated with (A) has the highest coherence value range from 0.6 to 1, whereas (B) and (D) explain the rate of damage according to the distance to the explosion. The slight damage is annotated with (D) showing the lowest interferometric coherence value 0.2. Maps Data: Google, image © 2021 Maxar Technologies.

Close modal
1

See https://disc.gsfc.nasa.gov/ (Last viewed ■).

3

See https://swp.ctbto.org/ (Last viewed ■).

1.
AIGA
(
2008
). “
Safe practices for the production of nitrous oxide from ammonium nitrate
,” Technical Report (AIGA); available at http://www.asiaiga.org/uploaded_docs/aiga%20080_16_safe_practices_for_the_production_of_nitrous_oxide_from_ammonium_nitrate.pdf.
2.
Archuleta
,
R. J.
,
Cranswick
,
E.
,
Mueller
,
C.
, and
Spudich
,
P.
(
1982
). “
Source parameters of the 1980 Mammoth Lakes, California, earthquake sequence
,”
J. Geophys. Res.: Solid Earth
87
(
B6
),
4595
4607
, .
3.
Attema
,
E.
,
Cafforio
,
C.
,
Gottwald
,
M.
,
Guccione
,
P.
,
Monti Guarnieri
,
A.
,
Rocca
,
F.
, and
Snoeij
,
P.
(
2010
). “
Flexible dynamic block adaptive quantization for sentinel-1 SAR missions
,”
IEEE Geosci. Remote Sens. Lett.
7
(
4
),
766
770
.
4.
BBC
(
2019
). “
China chemical blast death toll rises to 47
,” available at https://www.bbc.com/news/world-asia-47663023 (Last viewed 11/11/2020).
5.
BBC
(
2020
). “
Beirut explosion: Lebanon's government ‘to resign’ as death toll rises
,” available at https://www.bbc.com/news/world-middle-east-53720383 (Last viewed 17/10/2020).
6.
Blandford
,
R. R.
, and
Clauter
,
D. A.
(
1995
). “
Capability estimation of infrasound networks, AFTAC Report
,” Technical Report (AFTAC).
7.
Braun
,
A.
,
Fakhri
,
F.
, and
Hochschild
,
V.
(
2019
). “
Refugee camp monitoring and environmental change assessment of Kutupalong, Bangladesh, based on radar imagery of Sentinel-1 and ALOS-2
,”
Remote Sens.
11
(
17
),
2047
.
8.
Brown
,
D. J.
(
1999
). Summary of Infrasound Source Location Meeting: San Diego, November 9–10, 1998, Center for Monitoring Research Technical Report CMR-99/02, January.
9.
Ceranna
,
L.
,
Le Pichon
,
A.
,
Green
,
D. N.
, and
Mialle
,
P.
(
2009
). “
The Buncefield explosion: A benchmark for infrasound analysis across Central Europe
,”
Geophys. J. Int.
177
(
2
),
491
508
.
10.
Drob
,
D. P.
,
Emmert
,
J. T.
,
Meriwether
,
J. W.
,
Makela
,
J. J.
,
Doornbos
,
E.
,
Conde
,
M.
,
Hernandez
,
G.
,
Noto
,
J.
,
Zawdie
,
K. A.
,
McDonald
,
S. E.
,
Huba
,
J. D.
, and
Klenzing
,
J. H.
(
2015
). “
An update to the Horizontal Wind Model (HWM): The quiet time thermosphere
,”
Earth Space Sci.
2
(
7
),
301
319
.
11.
Drob
,
D. P.
,
Picone
,
J. M.
, and
Garcés
,
M.
(
2003
). “
Global morphology of infrasound propagation
,”
J. Geophys. Res.: Atmos.
108
(
21
),
1
12
, .
12.
Efthimiou
,
N.
,
Psomiadis
,
E.
, and
Panagos
,
P.
(
2020
). “
Fire severity and soil erosion susceptibility mapping using multi-temporal Earth Observation data: The case of Mati fatal wildfire in Eastern Attica, Greece
,”
Catena
187
(
July 2019
),
104320
.
13.
ElGabry
,
M.
,
Korrat
,
I.
,
Hussein
,
H.
, and
Hamama
,
I.
(
2017
). “
Infrasound detection of meteors
,”
NRIAG J. Astron. Geophys.
6
(
1
),
68
80
.
14.
Ford
,
S. R.
,
Rodgers
,
A. J.
,
Xu
,
H.
,
Templeton
,
D. C.
,
Harben
,
P.
,
Foxall
,
W.
, and
Reinke
,
R. E.
(
2014
). “
Partitioning of seismoacoustic energy and estimation of yield and height-of-burst/depth-of-burial for near-surface explosions
,”
Bull. Seismol. Soc. Am.
104
(
2
),
608
623
.
15.
Gaebler
,
P.
,
Ceranna
,
L.
,
Nooshiri
,
N.
,
Barth
,
A.
,
Cesca
,
S.
,
Frei
,
M.
,
Grünberg
,
I.
,
Hartmann
,
G.
,
Koch
,
K.
,
Pilger
,
C.
,
Ross
,
J. O.
, and
Dahm
,
T.
(
2019
). “
A multi-technology analysis of the 2017 North Korean nuclear test
,”
Solid Earth
10
(
1
),
59
78
.
16.
Garcés
,
M. A.
,
Hansen
,
R. A.
, and
Lindquist
,
K. G.
(
1998
). “
Traveltimes for infrasonic waves propagating in a stratified atmosphere
,”
Geophys. J. Int.
135
(
1
),
255
263
.
17.
Georges
,
T. M.
(
1972
). “
3D ray tracing for acoustic-gravity waves
,”
J. Acoust. Soc. Am.
51
(
1A
),
147
.
18.
GHD
(
2012
). “
Orica mining services: Report for Kooragang Island
,” March, available at https://apps.dtic.mil/sti/pdfs/ADA569462.pdf (Last viewed 24/12/2020).
19.
Gitterman
,
Y.
(
2010
). “
Sayarim Infrasound Calibration Explosion: Near-source and local observations and yield sstimation
,” in
2010 Monitoring Research Review: Ground-Based Nuclear Explosion Monitoring Technologies (Idc)
, pp.
708
719
, available at https://na22.nnsa.doe.gov/mrr/2010/PAPERS7.HTM (Last viewed 15/11/2020).
20.
Gitterman
,
Y.
,
Pinsky
,
V.
,
Amrat
,
A.-Q.
,
Darwish
,
J.
,
Mayyas
,
O.
,
Nakanishi
,
K.
, and
Hofstetter
,
R.
(
2005
). “
Source features, scaling, and location of calibration explosions in Israel and Jordan for Comprehensive Test Ban Treaty (CTBT) monitoring
,”
Israel J. Earth Sci.
54
(
4
),
199
217
.
21.
Goldstein
,
P. A. U. L.
, and
Snoke
,
A.
(
2005
). SAC availability for the IRIS community. Incorporated Research Institutions for Seismology Newsletter, 7(UCRL-JRNL-211140).
21.
Green
,
D. N.
, and
Bowers
,
D.
(
2010
). “
Estimating the detection capability of the International Monitoring System infrasound network
,”
J. Geophys. Res.: Atmos.
115
,
D18116
, .
22.
Green
,
D. N.
,
Vergoz
,
J.
,
Gibson
,
R.
,
Le Pichon
,
A.
, and
Ceranna
,
L.
(
2011
). “
Infrasound radiated by the Gerdec and Chelopechene explosions: Propagation along unexpected paths
,”
Geophys. J. Int.
185
(
2
),
890
910
.
23.
Groves
,
G. V.
(
1955
). “
Geometrical theory of sound propagation in the atmosphere
,”
J. Atmos. Terr. Phys.
7
,
113
127
.
24.
Hamama
,
I.
, and
Yamamoto
,
M.-y.
(
2021
). “
Infrasonic Earthquake Detectability Investigated in Southern Part of Japan, 2019
,”
Sensors
21
(
3
),
894
.
25.
Hersbach
,
H.
,
Bell
,
B.
,
Berrisford
,
P.
,
Biavati
,
G.
,
Horányi
,
A.
,
Muñoz Sabater
,
J.
,
Nicolas
,
J.
,
Peubey
,
C.
,
Radu
,
R.
,
Rozum
,
I.
,
Schepers
,
D.
,
Simmons
,
A.
,
Soci
,
C.
,
Dee
,
D.
, and
Thépaut
,
J.-N.
(
2018
). “
ERA5 hourly data on pressure levels from 1979 to present
,” available at https://www.copernicus.eu/en/access-data/copernicus-services-catalogue/era5-hourly-data-single-levels-1979-presen (Last viewed 10/12/2020).
26.
Hubbard
,
B.
,
Abi-Habib
,
M.
,
El-Naggar
M.
,
McCann
,
A.
,
Singhvi
,
A.
,
Glanz
,
J.
, and
White
,
J.
(
2020
). “
How a massive bomb came together in Beirut's port
,” The New York Times, available at https://www.nytimes.com/interactive/2020/09/09/world/middleeast/beirut-explosion.html (Last viewed 25/9/2020).
27.
Hutton
,
L. K.
, and
Boore
,
D. M.
(
1987
). “
The ML scale in Southern California
,”
Bull. Seismol. Soc. Am.
77
(
6
),
2074
2094
.
28.
Johnson
,
J. B.
(
2003
). “
Generation and propagation of infrasonic airwaves from volcanic explosions
,”
J. Volcanol. Geotherm. Res.
121
(
1-2
),
1
14
.
29.
Lingevitch
,
J. F.
,
Collins
,
M. D.
,
Dacol
,
D. K.
,
Drob
,
D. P.
,
Rogers
,
J. C. W.
, and
Siegmann
,
W. L.
(
2002
). “
A wide angle and high Mach number parabolic equation
,”
J. Acoust. Soc. Am.
111
(
2
),
729
734
.
30.
Lingevitch
,
J. F.
,
Collins
,
M. D.
,
Dacol
,
D. K.
,
Drob
,
D. P.
,
Rogers
,
J. C. W.
, and
Siegmann
,
W. L.
(
2002b
). “
A wide angle and high Mach number parabolic equation
,”
J. Acoust. Soc. Am.
111
(
2
),
729
734
.
31.
Livescience
(
2020
). “
Beirut blast was one of the biggest non-nuclear explosions ever
,” available at https://www.livescience.com/beirut-explosion-biggest-non-nuclear.html (Last viewed 8/11/2020).
32.
Martinis
,
S.
,
Caspard
,
M.
,
Plank
,
S.
,
Clandillon
,
S.
, and
Haouet
,
S.
(
2017
). “
Mapping burn scars, fire severity and soil erosion susceptibility in Southern France using multisensoral satellite data
,” in
2017 IEEE International Geoscience and Remote Sensing Symposium (IGARSS)
, pp.
1099
1102
.
33.
Massonnet
,
D.
,
Rossi
,
M.
,
Carmona
,
C.
,
Adragna
,
F.
,
Peltzer
,
G.
,
Feigl
,
K.
, and
Rabaute
,
T.
(
1993
). “
The displacement field of the Landers earthquake mapped by radar interferometry
,”
Nature
364
(
6433
),
138
142
.
34.
Molan
,
Y. E.
, and
Lu
,
Z.
(
2020
). “
Can InSAR coherence and closure phase be used to estimate soil moisture changes?
,”
Remote Sens.
12
(
9
),
1511
.
35.
Picone
,
J. M.
,
Hedin
,
A. E.
,
Drob
,
D. P.
, and
Aikin
,
A. C.
(
2002
). “
NRLMSISE-00 empirical model of the atmosphere: Statistical comparisons and scientific issues
,”
J. Geophys. Res.: Space Phys.
107
(
A12
),
SIA 15-1
SIA 15-16
, .
36.
Pierce
,
A. D.
(
1967
). “
Guided infrasonic modes in a temperature- and wind-stratified atmosphere
,”
J. Acoust. Soc. Am.
41
(
3
),
597
611
.
37.
Pierce
,
A. D.
, and
Posey
,
J. W.
(
1971
). “
Theory of the excitation and propagation of Lamb's atmospheric edge mode from nuclear explosions
,”
Geophys. J. Int.
26
(
1-4
),
341
368
.
38.
Pilger
,
C.
,
Gaebler
,
P.
,
Ceranna
,
L.
,
Le Pichon
,
A.
,
Vergoz
,
J.
,
Perttu
,
A.
,
Tailpied
,
D.
, and
Taisne
,
B.
(
2019
). “
Infrasound and seismoacoustic signatures of the 28 September 2018 Sulawesi super-shear earthquake
,”
Nat. Hazards Earth Syst. Sci.
19
(
12
),
2811
2825
.
39.
Pilger
,
C.
,
Gaebler
,
P.
,
Hupe
,
P.
,
Kalia
,
A. C.
,
Schneider
,
F. M.
,
Steinberg
,
A.
,
Sudhaus
,
H.
, and
Ceranna
,
L.
(
2021
). “
Yield estimation of the 2020 Beirut explosion using open access waveform and remote sensing data
,”
Sci. Rep.
11
(
1
),
14144
.
40.
Ramzy
,
A.
, and
Hernández
,
J. C.
(
2019
). “
Explosion at China chemical plant kills 64; Employees detained,” available at
https://www.nytimes.com/2019/03/22/world/asia/china-chemical-plant-explosion.html (Last viewed 7/10/2021).
41.
ReVelle
,
D. O.
, and
Kulichkov
,
S. N.
(
1998
). “
On Lamb wave propagation from small surface explosions in the atmospheric boundary layer,” available at
https://www.osti.gov/biblio/353177 (Last viewed 30/10/2021).
42.
Richter
,
C. F.
(
1935
). “
An instrumental earthquake magnitude scale
,”
Bull. Seismol. Soc. Am.
25
(
1
),
1
32
.
43.
Saito
,
H.
,
Yamamoto
,
T.
,
Nakajima
,
K.
,
Kuramoto
,
K.
, and
Yamamoto
,
M.-y.
(
2021
). “
Identification of the infrasound signals emitted by explosive eruption of Mt. Shinmoedake by three-dimensional ray tracing
,”
J. Acoust. Soc. Am.
149
(
1
),
591
598
.
44.
Schwaiger
,
H. F.
,
Iezzi
,
A. M.
, and
Fee
,
D.
(
2019
). “
AVO-G2S: A modified, open-source ground-to-space atmospheric specification for infrasound modeling
,”
Comput. Geosci.
125
,
90
97
.
44.
Smets
,
P. S. M.
,
Assink
,
J. D.
,
Le Pichon
,
A.
, and
Evers
,
L. G.
(
2016
). “
ECMWF SSW forecast evaluation using infrasound
,”
J. Geophys. Res.: Atmos.
121
(
9
),
4637
4650
, .
45.
Spaans
,
K.
, and
Hooper
,
A.
(
2016
). “
InSAR processing for volcano monitoring and other near-real time applications
,”
J. Geophys. Res.: Solid Earth
121
(
4
),
2947
2960
, .
46.
Stevens
,
J. L.
,
Adams
,
D. A.
,
Baker
,
G. E.
,
Xu
,
H.
,
Murphy
,
J. R.
,
Divnov
,
I.
, and
Bourchik
,
V. N.
(
2006
). “
Infrasound modeling using Soviet explosion data and instrument design criteria from experiments and simulations
,” Technical Report (Defense Technical Information Center); available at https://apps.dtic.mil/sti/citations/ADA526695.
47.
Stevens
,
J. L.
,
Divnov
,
I. I.
,
Adams
,
D. A.
,
Murphy
,
J. R.
, and
Bourchik
,
V. N.
(
2002
). “
Constraints on infrasound scaling and attenuation relations from Soviet explosion data
,”
Pure Appl. Geophy.
159
(
5
),
1045
1062
.
48.
Stump
,
B. W.
,
Pearson
,
D. C.
, and
Reinke
,
R. E.
(
1999
). “
Source comparisons between nuclear and chemical explosions detonated at Rainier Mesa, Nevada Test Site
,”
Bull. Seismol. Soc. Am.
89
(
2
),
409
422
.
49.
Tamkuan
,
N.
, and
Nagai
,
M.
(
2017
). “
Fusion of multi-temporal interferometric coherence and optical image data for the 2016 Kumamoto earthquake damage assessment
,”
ISPRS Int. J. Geo-Inf.
6
(
7
),
188
.
50.
Thompson
,
R. J.
(
1972
). “
Ray theory for an inhomogeneous moving medium
,”
J. Acoust. Soc. Am.
51
(
5B
),
1675
1682
.
51.
Tzouvaras
,
M.
,
Danezis
,
C.
, and
Hadjimitsis
,
D. G.
(
2020
). “
Small scale landslide detection using Sentinel-1 interferometric SAR coherence
,”
Remote Sens.
12
(
10
),
1560
.
52.
Watanabe
,
M.
,
Thapa
,
R. B.
,
Ohsumi
,
T.
,
Fujiwara
,
H.
,
Yonezawa
,
C.
,
Tomii
,
N.
, and
Suzuki
,
S.
(
2016
). “
Detection of damaged urban areas using interferometric SAR coherence change with PALSAR-2
,”
Earth, Planets Space
68
(
1
),
131
.
53.
Waxler
,
R.
,
Assink
,
J.
, and
Velea
,
D.
(
2017a
). “
Modal expansions for infrasound propagation and their implications for ground-to-ground propagation
,”
J. Acoust. Soc. Am.
141
(
2
),
1290
1307
.
54.
Waxler
,
R. M.
,
Assink
,
J. D.
,
Hetzer
,
C.
, and
Velea
,
D.
(
2017b
). “
NCPAprop—A software package for infrasound propagation modeling
,”
J. Acoust. Soc. Am.
141
(
5
),
3627
3627
.
55.
Wempen
,
J. M.
(
2020
). “
Application of DInSAR for short period monitoring of initial subsidence due to longwall mining in the mountain west United States
,”
Int. J. Min. Sci. Technol.
30
(
1
),
33
37
.
56.
Whitaker
,
R. W.
(
1995
). “
Infrasonic monitoring,” available at
https://www.osti.gov/biblio/760050 (Last viewed 11/11/2020).
57.
Xu
,
R.
,
Liu
,
J.
, and
Xu
,
J.
(
2018
). “
Extraction of high-precision urban impervious surfaces from Sentinel-2 multispectral imagery via modified linear spectral mixture analysis
,”
18
(
9
),
2873
.
58.
Zhang
,
N.
,
Shen
,
S.-L.
,
Zhou
,
A.-N.
, and
Chen
,
J.
(
2019
). “
A brief report on the March 21, 2019 explosions at a chemical factory in Xiangshui, China
,”
Process Saf. Prog.
38
(
2
),
e12060
.