Ionospheric Disturbances Produced by the 2019 Caribbean Meteoroid

1. Introduction

Bolides impact the Earth daily, ranging from small objects ($d<0.5$ m) that disintegrate upon atmospheric entry to larger ($d>100$ m), rare events capable of causing global disturbances. If their diameter is greater than 1 km the impact can be considered a global catastrophe [1,2]. Some famous modern cases are the Tunguska event 1n 1908 [3] and the Chelyabinsk event [4] in 2013.

In addition, the General Civil Protection Law in Mexico, amended in June 2014, incorporated near-Earth objects (NEOs) into the list of natural hazards phenomena [5]. The legislation mandates that the National Center for Disaster Prevention (CENAPRED) establishes responsibilities and the Mexican Space Agency to establish national action protocols in case of meteorites. Meteorite impacts and NEOS are also included in the UNDRR "Hazard Definition and Classification Review" [6], denoting their importance for the international community.

As bolides passes through the Earth’s atmosphere at hypersonic velocities, and posterior breakup, may produces turbulence, wave processes, etc. Thus leading to displacements of air, which produces Acoustic Gravity Waves (AGWs). AGWs can reach ionospheric levels and be detected as Traveling Ionospheric Disturbances (TIDs). Such events have been recently recorded in certain databases, such as the Center for Near Earth Objects Studies (CNEOS), Jet Propulsion Laboratory (JPL) and the Geostationary Lightning Mapper (GLM) publicly available at https://cneos.jpl.nasa.gov/fireballs/ and https://neo-bolide.ndc.nasa.gov/#/, respectively.

We picked our interest in one certain event, that we nicknamed “the Caribbean Meteoroid.” This meteoroid is the most energetic which appears in both databases and thus is an interesting case for studying the potential TID that could be produced. To do so, we collected Receiver INdependent EXchange (RINEX) data from nearby GNSS stations from the UNAVCO network (also publicly available 1) and developed a method to detrend the resultant TEC series and estimate the TIDs propagation speed.

The Figure 1 depicts a bright flash of light generated by the impact of a small Near-Earth Object (NEO) upon entering Earth’s atmosphere on June 22, 2019, over the Caribbean region. The flash, captured in a sequence of images, highlights the precise moment of the meteoroid’s entry and its interaction with the upper atmosphere just prior to disintegration. The image reveals a radiant expansion of energy, corresponding to the heat and pressure release caused by the meteoroid as it traverses the upper layers of the atmosphere. This visual context is critical for understanding the ionospheric disturbances induced by this event, which are discussed throughout the study.

Our approach involved advanced signal processing techniques, such as the Savitzky-Golay filter [7], to enhance the signal quality of the TEC time series. This allowed us to isolate wave-like features indicative of TIDs. We analyzed the spatial and temporal patterns of the disturbances to ascertain their velocities and directions of propagation, which are crucial for understanding the dynamics of ionospheric responses to meteoroid impacts [8].

Furthermore, we correlated these disturbances with the estimated impact energy derived from GLM light curves to establish a relationship between the meteoroid’s energy and the magnitude of ionospheric disturbances observed Maruyama et al. [9]. By comparing our findings with background ionospheric conditions, we aimed to delineate the influence of this specific meteoroid event against natural variability [10].

Our study not only sheds light on the immediate effects of meteoroid impacts on the ionosphere but also sets the stage for future investigations into similar phenomena, highlighting the need for more robust methodologies in detecting and analyzing TIDs caused by various atmospheric disturbances [11]. This research contributes to a broader understanding of the role of meteoroids in modifying the ionospheric environment and paves the way for advancements in predictive models for space weather events [12].

2. Data Processing

2.1. Meteor Location and Trajectory

The entry of the Caribbean meteor into the atmosphere was identified and studied its behavior and impact in the ionospheric region. To do so we used data from the GLM [13] and the CNEOS of the JPL of the NASA. We used the interactive database of both projects, available at https://neo-bolide.ndc.nasa.gov/#/ and https://cneos.jpl.nasa.gov/fireballs/. From CNEOS database, we found basic physical parameters such as event coordinates, altitude, total velocity and its components, total radiated energy and calculated impact energy. These data are shown in Table 1. On the other hand, from GLM we got the meteor light curve, accurate position of the fragments and duration event. The use of both databases allowed us to cross-verify the detection data, providing complementary information about the meteor’s properties due to the independent nature of their measurements.

According to the CNEOS database records, the total energy is about 1.4% of the Chelyabinsk bolide, and fragmentation occurred above the open sea, at afternoon or near sunset depending on the closest mainland available: from Dominican republic, Puerto Rico, the lesser Antilles, Venezuela, Guyana, French Guyana and Suriname.

3. RINEX Data

In order to study the ionospheric behavior, we use the RINEX data for event from 38 stations that surrounded the meteor detection coordinates for the day of the event. Since the event occurred near UT midnight, we also obtained GNSS data for the next day in case of there would be TIDs detected several hours after the fragmentation. Using a software developed by Gopi K. Seemala [14], publicly available at https://seemala.blogspot.com/, we computed the slant TEC (sTEC) and vertical TEC (vTEC) for different GNSS satellites, each one identified with a PseudoRandom Noise code (PRN). The behavior of the TEC curve is due to many factors, including the Earth’s rotation, solar activity, etc. TID’s and wave-like curves are not as prominent and are difficult to identify. Due to the sample size, the detrending process necessitates manual intervention to ensure precise data interpretation [15]. In Appendix A1 shows the stations names, coordinates and proper citations.

3.1. Detrending Process

For detrending our data we used a method developed for detecting plasma bubbles in the equatorial region [16], but proved to be effective for detecting AGWs and TIDs. With this method we are able to infer the trend from our TEC data using a Savitzky-Golay filter [16,17]. However, the results of this filtering are sensitive to the parameters we explored or use for such filtering, thus altering significantly the quality of the detrended signal. The TEC time series are assumed to be an additive combination of a signal and a trend [15,18]. Thus we tested the method of Pradipta et al. [16] over the sum of a known signal and a known trend and find the parameters of the Savitzky-Golay filter that recover in the best way the original signal. We used a test signal with the purpose to find out the most appropriate parameters for the filter and we show its functional form in Appendix B.

3.2. Comparison with ARTU data for Chelyabinsk meteor

Our final test was to compare our resulting time series with the results of another, more documented event: the 2013 Chelyabinsk meteor impact. Yang et al. [4] obtained TEC time series by analyzing the coherence spectrum of TEC measurements with many GNSS stations distributed in Russia, United States, and Japan, with the ARTU station being the closest to the meteor impact location. Then, using the most prominent detected frequencies (4.0-7.8 mHz and 1.0-2.5 mHz), they reconstructed the TEC series for PRN 18 and 26. This methodology aligns with previous findings that emphasize the significance of frequency analysis in identifying ionospheric disturbances [19].

Thus, we obtained the RINEX data for the ARTU station and detrended the data using the method described in this paper for the ARTU-PRN18 and ARTU-PRN26 station-satellite links. The results are shown in Figure 2. In this case, we have noticed that the positions of the resulting AGWs align, so we can be confident in our method, which is consistent with the work of [4] with TIDs periods of around 20 min. However, it is also noteworthy that the shape of the detrended time series is not the same due to the fact that it is actually the superposition of all the frequencies involved [20].

Additionally, we already knew where the produced ionospheric perturbation occurred. Finding and distinguishing ionospheric perturbations using this method can be a difficult task, but it has proven to be an accessible confirmation method. This aligns with the conclusions of previous research that highlights the utility of detrending techniques in isolating significant disturbances from background noise [21].

4. Meteor Physical Properties

Using the estimation of the meteor energy and velocity from Table 1, the trajectory angle estimated in Appendix A from the velocity components of ${25.6}^{\circ }$ and assuming the meteor has a density of 4000 kg m⁻³, roughly consistent with the composition of chondrites [22], we estimated the initial velocity, initial mass and entry angle of the Caribbean meteor, enlisted in Table 2, by using a system of equations that models the dynamics of a meteoroid passing through Earth’s atmosphere [23]. Solving numerically this set of equations we modeled the Caribbean meteor evolution as passing through Earth’s atmosphere. In Figure 3, Figure 4 and Figure 5, we show the evolution of meteor parameters as it passes through the atmosphere.

In Figure 3, we see that after breakup mass and velocity drops fast to zero, this is due that after breakup these graphs show the mass and velocity of the fragment cloud, which is more susceptible to dissipation, explaining the fast mass drop, each individual fragment into the cloud is easily slowed down by the atmosphere, explaining the velocity drop. The model lost trace of the fragment cloud mass slightly below 30 km. On the other hand, the trajectory angle changes with height linearly in about one degree similar to the suggested by [20,24,25].

The meteor diameter and density are shown in Figure 4, we see for density in left panel that it drops drastically at 40 km approximately due that the density considered after breakup is for the fragmentation cloud, where all the meteor fragments are scattered in a greater volume, for the same reason in the right panel the diameter increases from the same height [26].

Finally, we obtained the meteor energy and pressure, shown in Figure 5. We observe the meteor energy as a function of height, which decreases after breakup and the burst, since most of this energy was released after the burst, which is behavior similar to that suggested by [3,27]. The deposited energy is shown in top panel of Figure 5, where it is discernible a similar behavior however is similar but it has maximum around a height of 40 km. Finally the pressure is shown in right panel, where the meteor stagnation pressure is also displayed. When the meteor pressure exceeds the stagnation pressure it increases dramatically until a maximum at $h\approx 25$ km.

5. Estimated meteor trajectory

Using the data obtained from the model explained in the former section, we estimated the meteor trajectory. The data extracted from the model is shown in Table 3.

With these data, the meteor trajectory can be estimated as follows:

$$\begin{array}{cc}\hfill L-L_0& =msin\alpha {\int }_{h_0}^h{\displaystyle \frac{vcos\theta }{R_E+h}}dh,\end{array}$$

(1)

$$\begin{array}{cc}\hfill \lambda -{\lambda }_0& =mcos\alpha {\int }_{h_0}^h{\displaystyle \frac{vcos\theta }{R_E+h}}dh.\end{array}$$

(2)

Where L and $\lambda$ correspond to the latitude and longitude as a function of height h, respectively, the 0 sub-index correspond to the position where measurements started, in this case $h_0=70$ km, v is the meteor velocity, $\theta$ is the trajectory angle, $R_E$ is the radius of Earth and $\alpha$ is the trajectory angle respect to the equator, assumed to be constant since Coriolis force is negligible for this case, and computed from the components of the meteor velocity in Table 1. Following the results described in §4, we found that a linear relation between time and height, being m the proportionality constant. For heights above 70 km until the Ionospheric Piercing Point (IPP), at 350 km we assumed that the meteor was moving at constant speed equal to the terminal velocity and we supposed a linear relation between trajectory angle between height and trajectory angle $\theta$.

6. Ionospheric Background and sTEC Time Series

Ionospheric perturbations can take place due to space weather conditions, with an origin from solar events affecting the Earth environment by radiation or charged particles. So, in order to discard such events we investigated the space weather conditions in the day each event occurred. We investigated the solar wind parameters in the event day and the previous 7 days, the x-ray flux and the Dst index at the event date.

6.1. Dst Index

Measurements of the Dst Index were obtained from WDC for Geomagnetism, Kyoto DST index service (https://wdc.kugi.kyoto-u.ac.jp/dstdir/). From these data, we found out that Dst index remains practically constant around zero, which confirms that there were not geomagnetic storms that could prevent us from detecting TIDs in our GNSS data, see Figure 6.

6.2. Solar Wind

The solar wind is part of the space weather events that we should take in account (in order to confirm if TIDs produced by the meteor passage are viable to be detected or no) is unusual solar wind behavior or even coronal mass ejections. To do so, we collected data from the Deep Space Climate Observatory (DSCOVR) of solar wind speed, temperature, density, dynamic pressure and solar magnetic field, in order to monitor the solar wind behavior in the day of the event and the previous 6 or 7 days.

Potential shock waves produced by rapid wind transitions were considered reaching slow wind, which also may produce TIDs and thus contaminating our GNSS data, and may not be distinguishable from TIDs produced by the meteor passage. To check the presence of such shock waves, we estimate the relative velocity of solar wind of one day respect the previous one. The results are shown in Figure 7, where we observe that this velocity is almost constant and gives no chance for a fast stream reaching a slow that could produce a shock wave.

6.3. X-Rays Flux

X rays can be produced in solar flares impacting directly the ionospheric layer. A sudden increase of the flux or sudden variations may be a source of AGW or giving appearance of AGW. We investigated the X ray flux at the event date and the previous week, and collected data from the NOAA Space Prediction Center. A thorough review confirmed the absence of type C, M, or X solar flares capable of impacting Earth’s ionosphere during the period in question.

6.4. Solar Terminator

When a meteoroid enters the atmosphere and experiences subsequent breakup or burst near the solar terminator, the influence of the terminator’s passage on ionospheric conditions warrants careful consideration. Acoustic gravity waves (AGW), observed near the terminator region, are known to induce significant disturbances in the ionosphere [18]. Therefore, disentangling meteoroid-induced effects from terminator-related phenomena in such scenarios requires a comprehensive analysis that accounts for the potential contributions of both sources [28]. Figure 8 shows that meteor burst occurred outside the solar terminator, but some GNSS stations (KOUG and KOUR) are located close to (20 minutes after the fragmentation). In an effort to distinguish between the effects of the meteor and the solar terminator we detrended GNSS data for the mentioned stations at dates prior of the meteor fall, those dates are close enough to guarantee that the sunset occurs at different times. According to the previous, we can isolate the effects of the solar terminator and analyze its effects to check how they could affect time series. Some results are shown in Figure 9 for station KOUG, PRNs 2 and 12. In these cases we see that the solar terminator perturbs the ionosphere, and its effects can bypass the detrending process, since they consist in wave-like variations of TEC. However, these perturbations are typically weak and can be distinguished from those generated by meteor passage, allowing us to confidently conclude that they do not significantly affect our detections. Even more, comparing time series from previous days it can be helpful us to discriminate between time series with detected TID from others where only noise is recorded since some patterns in the time series can be detected. The time series that are similar with the previous days can be classified as non detections.

6.5. sTEC Time Series

We processed the satellite-receiver time series collected from all stations, analyzing each series separately to identify wave-like features through visual inspection. Caution is required when drawing conclusions regarding the presence of TIDs or AGWs associated with the bolide, as ionospheric perturbations induced by the meteor passage may be confounded with other types of ionospheric disturbances. A significant sample of about 10% of all the time series are depicted in Figure 10 and Figure 11.

To facilitate the organization of data and trace spatially the propagation of TIDs, we divided the GNSS stations into 5 sets, roughly based in their location. The spatial distribution of sets and stations is shown in Figure 12.

After analyzing data related to TIDs, the analysis indicates that the stations detecting TIDs are part of sets 3, 4 and 5, which are located southeast of the meteor detection site. It is suspected that TIDs were generated at the position where the meteor crossed through the IPP and propagated away from the source. To estimate such position, we reconstructed the meteor trajectory using the results of the flight equations [24], but these results only were useful up to a height of 70 km. To reconstruct the rest of the trajectory, we observed that before breakup the meteor was moving at a terminal velocity of $20.6$ kms⁻¹, and then extrapolated the linear correlation between trajectory angle and height (see, for example, Figure 3). The resultant trajectory is also shown in Figure 12, as well as the positions of breakup and the meteor crossing the IPP (blue star), from where the TIDs should have been generated (TIDs source from this point). Stations from set 3 and a few from set 4 were close enough to detect TIDs as they passed through. In fact, TIDs were detected for stations located at the south of the meteor estimated trajectory. Set 1 stations were the farthest of the TIDs source. If stations of this set detected anything, the TIDs amplitude after traveling such distances should have been diminished, and become indistinguishable from local noises. The same applies for stations from set 2, but is evident that we did not find TIDs from this set of stations, since the distance to the TIDs source is much shorter, and is also evident that the stations from set 4 and 5 detected TIDs farther from the source. To explain this situation, an independent detection method should be used as corroboration. In this context, we used ROTI index for high resolution data in order to study this.

6.6. High Resolution Data and ROTI Estimation

We employ high resolution GNSS data (1 Hz) from a subset of stations enlisted in Table A1 of Appendix A. From these data, we recalculated dTEC in order to see if the dTEC changes correlated with the meteoroid event. By employing the ROTI method, we were able to quantify the rapid fluctuations in TEC that typically accompany TIDs. ROTI is a useful tool for assessing ionospheric irregularities, as it highlights variations that may not be evident in standard TEC measurements [29].

The dTEC time series were analyzed before, during, and after the impact to capture the temporal evolution of the disturbances. This analysis allowed us to observe significant spikes in dTEC that corresponded with the expected arrival time of the TIDs generated by the meteoroid’s atmospheric entry. The use of high-resolution data enhanced our ability to detect these transient events, providing a clearer picture of their propagation characteristics [30,31].

Furthermore, by applying statistical methods to the ROTI values, thresholds were established to identify significant disturbances linked to meteoroid impacts. This approach not only improves our detection capabilities but also contributes to understanding the broader implications of such events on ionospheric dynamics [31]. Our findings emphasize the importance of using high-resolution GNSS data and advanced analysis techniques to accurately capture the intricate behavior of the ionosphere during meteoroid-induced disturbances.

Nevertheless, with such data, ROTI may be computed. It is defined as the standard deviation of the rate of change of TEC (ROT) within an interval of 1 minute and is computed forward in time through moving average. Low resolution data provides insufficient information to compute satisfactorily the ROTI, but high resolution data solves this issue. In Figure 13, we compare the computed ROTI for station CN40 for satellites with PRN 2, 5, 6, 9, 12 and 15 (which seems to be the most probable TIDs detections) against the detrended data.

Scintillation indices such as S4 and ${\sigma }_{\phi }$ are not available in this study, but ROTI and Multipath (MP) are correlated with scintillation under certain circumstances that are fulfilled in this study [32]. MP is computed with the following the equation:

$$MP={\displaystyle \frac{1}{N}}\sum _{i=1}^N{\left(TE{C}_i-\overline{TEC}\right)}^2,$$

(3)

where N is the number of data points, $TE{C}_i$ is the total electron content at time i, and $\overline{TEC}$ is the mean TEC over the observation period. This approach helps to quantify the fluctuations in TEC that arise from multipath effects, which can occur when satellite signals reflect off nearby surfaces before reaching the receiver.

In our analysis, we focused on periods of high ionospheric activity associated with the meteoroid impact. By correlating MP with ROTI, we aimed to enhance our understanding of how these metrics interact and their implications for ionospheric disturbances. Previous studies have shown that higher multipath effects often coincide with increased scintillation activity, particularly during geomagnetic storms and other significant ionospheric events [33,34].

These insights are crucial for interpreting the overall impact of the meteoroid on the ionosphere and for refining models that predict ionospheric behavior under similar circumstances. Our findings suggest that monitoring both ROTI and MP can provide a more comprehensive view of the ionospheric response to transient disturbances. However, on this occasion, we did not observe any evidence of perturbations in the ROTI index.

Finally, dTEC and the W index were also calculated for each station on the event day. The W index, indicative of ionospheric perturbation levels, did not exhibit any significant fluctuations, aligning with the stable behavior observed in the Dst index.

7. Discussion and Conclusions

GNSS data from a subset of stations within the UNAVCO network was collected, focusing on the period surrounding the detection of the Caribbean bolide event. Data acquisition emphasized the day of the event and the subsequent day. Additionally, GLM data, specifically total energy estimates derived from light curves, provided measurements enabling the determination of the bolide’s position and the estimation of its trajectory at the point of fragmentation.

For the GNSS data, we adapted the method from Pradipta et al. [16] to detrend the resultant TEC curves, effectively removing the influences of Earth’s rotation and solar activity (including diurnal and seasonal variations) to make the wave-like features more pronounced. We estimated initial parameters of the meteor, such as velocity, mass, and trajectory angle, assuming a density of 4000 kg m⁻³, which is comparable to that of chondrites. We also analyzed how these parameters—velocity, mass, density, diameter, and trajectory angle—changed as the meteor descended through the atmosphere, allowing us to derive energy and pressure. Notably, we found no evidence of intense solar activity on the day of the meteor’s fall or in the preceding days. Indicators such as unusual fluctuations in x-ray flux, spikes in solar wind speed, and significant variations in the Dst and W indices remained within normal levels. Only two (or perhaps three) stations were located near the solar terminator at the time of meteor detection, differing by about an hour. To isolate the solar terminator’s effects at these stations, we adjusted time series from the days leading up to the event and discovered that while solar terminator effects were present, they were weaker and had minimal impact on our detections.

The findings of this study confirm that the 2019 Caribbean meteoroid event generated detectable ionospheric disturbances, evident in the detrended sTEC signal. This phenomenon, consistent with observations from events such as the 2013 Chelyabinsk meteor airburst, underscores the capacity of extraterrestrial objects to induce significant perturbations in the upper atmosphere. Comparisons with the Chelyabinsk event reveal similarities in the propagation characteristics of ionospheric disturbances and the velocity of Acoustic Gravity Waves. However, the Caribbean event exhibited lower overall energy.

A key contribution of this work is the validation of adjustment methods for analyzing meteoroid-induced ionospheric disturbances. The study demonstrates the efficacy of employed tools, such as the Savitzky-Golay filter, in identifying Acoustic Gravity Waves and Traveling Ionospheric Disturbances, consistent with the work of [4] with similar periods of around 20 min for these AGWs. Furthermore, data integration from networks like UNAVCO and databases such as CNEOS and GLM facilitated the development of a more accurate model of the ionospheric impact. However, challenges remain in detecting and analyzing these phenomena due to data sensitivity limitations and inherent variations in Total Electron Content.

Having ruled out the presence of extraneous Traveling Ionospheric Disturbances and Acoustic Gravity Wave phenomena, the observed wave-like features within the time series are confidently attributed to ionospheric perturbations resulting from the Caribbean meteor’s atmospheric passage. In summary, this study confirms the significant role of meteoroids in generating ionospheric disturbances, validating the effectiveness of advanced monitoring tools such as the Savitzky-Golay filter, and informing the development of enhanced detection methods for space weather research. Given the potential atmospheric and security implications of space object impacts, continued research in this domain is essential to improve our understanding and response capabilities for future events. Future publications will encompass estimations of the ionospherically perturbed area, analysis of TID velocity propagation, and wavelet analysis.

The insights derived from this analysis can be used to refine predictive models of ionospheric behavior, particularly in response to transient space objects, thereby enhancing space weather forecasting capabilities. This study establishes new avenues for investigating the complex interactions between meteoroids and the Earth’s atmosphere. Future research directions include analyzing similar events across varying latitudes and geophysical conditions, as well as refining models of Acoustic Gravity Wave propagation through different atmospheric layers. The implementation of more robust methods for data filtering and TID detection, such as incorporating the Rate of TEC Index (ROTI) to identify ionospheric disturbances with greater resolution, is also advised. Given the multifaceted implications of meteoroid impacts on ionospheric dynamics, encompassing atmospheric and security-related concerns, continued advancements in detection methodologies and predictive modeling efforts within this domain are imperative.