Bolide Energetics and Infrasound Propagation: Exploring the 18 December 2018 Bering Sea Event to Identify Limitations of Empirical and Numerical Models Open Access

Infrasound observations provide an important forensic clue for understanding meteoroids, particularly their energetics (Silber et al., 2018). The most widely used method for estimating bolide energetics from infrasound observations is based on a regression model that relates the dominant signal period to the equivalent nuclear yield, as determined from observations of atmospheric nuclear tests (Revelle, 1997). The popularity of this model is due, in part, to the fact that infrasound period is relatively unperturbed by atmospheric propagation effects in comparison with amplitude, the other indicator of source yield. However, it is recognized that there are important source differences between bolides and atmospheric nuclear tests. In particular, two key distinctions are that bolides have substantial mass and are not strictly point sources, although they can appear and be treated as such at sufficiently long range. Despite these discrepancies, the infrasound period–yield model has been applied to assess the energetics of a number of bolides (e.g., Revelle, 1997; Ens et al., 2012; Pilger et al., 2020).

Alternatives to using regression models require simulating infrasound propagation over long ranges, a problem that is extremely challenging because such simulations require specifications of the state of the atmosphere between the source and receiver at high spatial and temporal resolution from the ground to the upper atmosphere (Drob et al., 2015). A large body of recent work has shown that infrasound propagation is sensitive to atmospheric phenomena with spatial and temporal scales that are unresolved in state‐of‐the‐art specifications, in particular gravity waves (e.g., Hedlin and Drob, 2014). At distances of less than 200 km, the effect of gravity waves has been used to explain stratospheric observations within shadow zones (e.g., Hedlin and Drob, 2014). It is not surprising, therefore, that infrasound propagation models often fail to predict observations. However, observations that highlight systematic discrepancies between model simulations and observations, or that capture a key discrepancy at unprecedented resolution or scale, are important exemplars for motivating the development of numerical propagation models. In this article, we highlight such a discrepancy between predicted and observed transmission loss along an extremely long range (∼8000 km) that is sampled by two dense infrasound networks.

This study is based on infrasound observations of a large bolide that occurred over the Bering Sea on 18 December 2018, and was detected and located by satellite observations (Borovička et al., 2020). The event was reported on the NASA fireball website (see Data and Resources) with an origin time of 23:48:20Z, and location at 56.9° N and 172.4° E with altitude of peak brightness at 25.6 km. Using regression models between optical luminescence and yield NASA reported an impact energy of 173 kT Trinitrotulene (TNT) equivalent (1 kt of $TNT=4.184×1012J$⁠). This estimate would make this event the second largest recorded on the International Monitoring System (IMS) infrasound array, following the Chelyabinsk event (e.g., LePichon et al., 2013). The infrasound observations from the Bering Sea event on the IMS network have been reported by Pilger et al. (2020), who estimate a yield of 50 kT using an empirical period–yield law, which we discuss more in the following. In their paper, Pilger et al. (2020) note the unusual discrepancy with the yield derived from optical measurements, as compared with nine other events. This article is also based on infrasound signals from this event, but we incorporate additional high‐density regional observations. In addition, we note that this study is not focused on understanding the event itself but on exploiting the unique set of observations to assess empirical and numerical models, and in particular to highlight and quantify a key discrepancy between observations and models that it is not currently possible to resolve.

The Bering Sea bolide was observed on both IMS arrays and high‐density infrasound networks (Figs. 1, 2, 3). Although the IMS observations span a wide range of azimuths and distances, there are a particularly large number of observations toward the northeast due to the presence of the Transportable Array (TA; Busby and Aderhold, 2020) and Central and Eastern United States Network (CEUSN; Hedlin et al., 2018). This high density of observations over a long range makes this a unique set of observations.

Analyst measurements of the infrasound observations at the IMS arrays are provided in Table S1. These measurements, and the data shown in Figure 2, highlight that signals toward the west generally have lower signal‐to‐noise ratios (SNRs) and higher dominant frequencies than signals to the north and east. The distribution of SNR’s is consistent with the seasonal pattern in stratospheric winds, which drive infrasound detection capability at global scales and are favorable for eastward propagation in the northern hemisphere wintertime (Negraru and Golden, 2017; Green et al., 2018). Using the optical location and origin time as a type of ground truth (Green et al., 2010), we compute celerities (we follow the infrasound convention of using the term “celerity” to refer to the horizontal range divided by the travel time) for all arrivals (Table S1) and find that signals propagating to the east show fast celerities above 300 m/s, whereas signals observed west of the source always have lower celerities than those to the east. We observe significant variation in the observed periods (Table S1), ranging from 2.4 to 21 s. An analysis of the data has found that this variability is driven by site noise conditions, which corrupt the measurement in which the site noise is high or signal strength is low. This is illustrated in Figure 2 for the observations at the two closest IMS stations, I44RU and I53US, located at distances of 1025 and 2278 km west and east from the source. I44RU shows no observable signal above the noise for frequencies below 0.3 Hz. In contrast, I53US, even if it is further than I44RU, has a high SNR, and the signal is above the noise almost in the entire bandwidth, with the peak power around a frequency of 0.07. The I53US noise spectrum also exhibits a prominent microbarom peak—usually a good indicator of low background noise conditions—whereas I44RU does not exhibit such a peak. For consistency purposes we determined the dominant periods from the zero crossing of the 0.02–1 Hz filtered beam for the signals with high SNR. We also provide in Table S1 the dominant periods of the observed signals for the other arrays, which were determined in a different bandwidth (usually 0.2–1 Hz), but those periods are related to the site noise conditions rather than the bolide energetics. Several other explanations have been proposed to explain variations in the observed signal periods from bolides, such as departures from an explosive spherical source (Silber et al., 2009), Doppler shifts caused by winds along the propagation path (Ens et al., 2012), dispersion effects due to propagation at different ranges (Revelle, 1976), height of burst effects (Herrin et al., 2008), and site noise that can mask the true spectrum of the signals (Bowman et al., 2005). In the IMS dataset analyzed here, it is clear that at least some variation is due to the local noise conditions that affects the bandwidth of the observations.

The observations on the TA and CEUSN networks are distributed over a swath of azimuths from 30° to 75° relative to the location published by NASA (Fig. 1). We made measurements of infrasound signal amplitude on a subset of 79 stations with high SNR (Table S1) and determined that, similarly to the IMS network, the SNR is driven by site noise. The subset of stations, shown in Figure 3, exhibit clear infrasound arrivals with consistent celerities over a long propagation path between 0.29 and 0.315 km/s. Arrival times of signals on TA stations in Alaska are the best fit by a consistent celerity of 0.313 km/s, whereas arrival times of signals on CEUSN stations in the eastern United States are the best fit by a celerity of 0.305 km/s (Fig. 3). The consistent move out of the signals strongly suggests propagation in a persistent duct across Alaska, with a slight change in propagation occurring between 4000 and 6000 km due to slightly lower celerities for arrivals on CEUSN stations. This set of observations provides an unusual opportunity to explore the propagation of infrasound over a long propagation path (out to 8500 km) with high‐density observations over much of the path except between 3500 and 6000 km.

Assuming a log normal distribution of the observed periods at the 8 sites (out of a total of 13) with good SNR, a mean period of 14.8 ± 3.74 s in formula 1 gives a range of 2–84 kt with a mean of 21.4 kt. Because blast‐wave measurements are dependent on the surrounding medium, some authors correct for the lower air density at the altitude of the blast (Herrin et al., 2008; Antolik et al., 2014). The correction used is W exp(−Z/H), in which W is the yield, Z is the altitude (in km) of the explosion, and H is the scale height of the atmosphere, taken here as 8 km. Such a correction would lower the yield of the event by more than one order of magnitude to about 0.5 kt. It is worth noting that most of the authors who used the period–yield relationship to infer bolide energy (Silber et al., 2009; Pilger et al., 2020) do not apply this additional correction. Regardless of how the factor of two is applied, or whether additional corrections are made, it is clear that all the estimates from the signal periods are significantly lower than the NASA published yield of 173 kt. If we assume that the period–yield equation is correct and applicable to bolide energy, the higher energy estimate would require a period of 31 s (neglecting the altitude correction), which is not observed at any of the global sites, regardless of the observed SNR and bandwidth of the signals.

The coefficients A, B, and C are shown in Table S2, and the units for P, W, and R are Pascals, kilotons, and km, respectively, with the exception of Clauter and Blandford (1998), which uses degrees. The term controlling the rate of falloff with distance, C, implicitly includes geometric spreading and atmospheric absorption effects. For pure geometrical spreading, useful reference values of C are 1 (for spherical spreading) and 0.5 (for cylindrical spreading). In the published literature, values of this term range between 0.5 and 1.76 (Table S2). As shown in Figure 4, none of the relationships fits the dataset in its entirety, though there is a general agreement between the observations and the formulas. Using the NASA optical yield estimate, the Russian downwind formula fits the TA data up to around 2500 km; the Clauter and Blandford (1998) formula fits some of the CEUSN data, but even the high exponent of 1.76 in Blanc et al. (1997) is not sufficient to fit observations adequately. A linear regression between logP and log R shows that the best fits for the TA dataset only is obtained for an exponent of 2.2, whereas the best fit for the CEUSN dataset is obtained for an exponent of 2.51. If we use the entire dataset in the regression, we obtain an exponent of 2.42 (blue line in Fig. 5). This is much higher than the exponents in the earlier equations. Clearly, for this event, the application of these empirical relations to estimate yield would be problematic, because the distance exponent underpredicts the observed attenuation with range.

The results shown so far highlight that the observed celerities of the arrivals on TA and CEUSN stations suggest the presence of a persistent stratospheric duct, with stable properties out to at least 3500 km, and some variation after that. Furthermore, the amplitudes falloff at a rate much higher than expected for geometric spreading in a waveguide. Given that the absorption of sound in the stratosphere is low (Sutherland and Bass, 2004), the falloff is particularly surprising. To model the infrasound propagation, we extract profiles from the Ground‐to‐Space model along the 60° transect shown in Figure 1. The profiles suggest the presence of a strong stratospheric duct that gradually weakens at a range of around 3500 km and then reappears after 6000 km.

We use a parabolic equation (PE) model (Waxler et al., 2017), at a frequency of 0.1 Hz (corresponding to the average period measurement in Fig. 2), with 4 Pade coefficients, and including atmospheric absorption, to model the propagation along the 60° transect using a range‐dependent atmospheric model. The results, shown in Figure 5, indicate that the PE model predicts propagation in a stratospheric duct out to a range of ∼3500 km (covering the range of distances of TA stations), changing to a duct in the mesosphere and thermosphere from 3500–6000 km, and returning to stratospheric ducting at distances >6000 km (where the CEUSN stations are located). The estimated transmission losses, including absorption, are less than transmission losses predicted by pure spherical spreading with no absorption. It is important to emphasize that the PE simulations do not account for a spherical earth model, which has been shown to have a significant effect on predictions of geometric spreading in seismology (e.g., Yang et al., 2007). Furthermore, the PE model assumes a flat ground, thus does not account for scattering and loss of energy associated with each ground reflection.

To obtain further insight on propagation at global distances for a source at 25.6 km altitude, we simulate raypaths using raytracing equations derived for spherical coordinates and accounting for a 3D atmosphere (Blom, 2019). The results (shown in Fig. S1) indicate that the observed celerities are consistent with stratospheric returns, and that thermospheric returns are associated with celerities significantly below the observed values. The results further highlight the change in propagation characteristics at distances beyond 4000 km, with no raypaths predicted to hit the ground at larger distances. In this case, it appears that full‐wave models are necessary to predict the signals observed at the CEUSN station. A key conclusion of the full‐wave propagation modeling is that it does not predict the observed amplitude decay with range that decays faster, but consistently predicts a decay between pure spherical and cylindrical spreading.

Our analysis of this dataset illustrates the challenges in estimating a consistent energy release for a bolide event using both empirical and propagation‐based methods. All the methods presented here, empirical period–yield, empirical pressure–yield, and modeling‐based approaches, show limitations in their applicability.

The main challenge of using the period–yield relationship is obtaining a consistent period estimate. Although dispersion effects due to propagation at different ranges are expected, these effects are thought to be small (Revelle, 1997). Our analysis suggests that site noise is the dominant effect, and accurate determination of the period can be dependent on the SNR as a function of frequency. Agreement between the NASA optical estimate and the infrasound data would require a mean dominant period of 31 s (calculated by solving for the period in formula 1 using a yield of 173 kt), which is not observed at any of the sites.

The signal amplitudes at the TA and CEUSN provided us with a unique opportunity to examine the attenuation using observations across a large range of distances. The observed attenuation exceeds the highest published attenuation coefficient of 1.76 (Table S2), with the best fitting attenuation coefficients of 2.21–2.51 observed here. Such coefficients are larger than the previously published relationships and would result in a lower yield. Because of the discrepancy between our observations and published models, we do not attempt to estimate a yield from these amplitudes.

The observed change in attenuation coefficient from 2.21 to 2.51 from the TA to the CEUSN is consistent with propagation modeling, which shows changes at distances beyond 3500 km, around the maximum range of the TA observations. Previous work by Blom et al. (2018) has developed stochastic methods to account for transmission losses to infer source yields. The results presented here suggest that further empirical assessment of transmission loss with range is necessary to understand the limitations of models. In particular, this article highlights a unique set of observations over a long range, in which the amplitude attenuation does not match the predicted transmission losses and exceeds previous empirical attenuation rates. Because the attenuation is higher than predicted, possible discrepancies may relate to the effects of losses due to scattering off topography, losses due to out‐of‐plane effects, or the effects of scattering due to fine‐scale atmospheric structure that is not resolved in state‐of‐the‐art atmospheric models. Models capable of accounting for these effects are not yet computationally feasible at these scales. However, the observations in this article provide a useful benchmark for future propagation codes and atmospheric models.

This article uses data from the United States Transportable Array experiment and the Central and Eastern United States Network, which can be obtained from the Incorporated Research Institutions for Seismology (IRIS) Data Management Center at www.iris.edu. In addition, we use data from the International Monitoring System (IMS) of the Comprehensive Test Ban Treaty Organization. IMS data from the U.S. stations can also be obtained from IRIS Data Management Center.

Some plots in this article were made using the Generic Mapping Tools 4.5.6 (www.soest.hawaii.edu/gmt; Wessel and Smith, 1998). The 18 December 2018 Bering Sea event was reported on the NASA fireball website (https://cneos.jpl.nasa.gov/fireballs). All websites were last accessed in December 2021. Supplemental material includes Figure S1, showing raypaths calculated along a 60° transect. Supplemental Table S1 contains measurements made on the IMS arrays and TA stations. Supplemental Table S2 contains the coefficients used for the pressure formulas.

The authors acknowledge that there are no conflicts of interest recorded.

The authors would like to acknowledge the peer review and helpful comments from Dr. Mark Woods.