# A Numerical Approach to Study Ablation of Large Bolides: Application to Chelyabinsk

Authors: Josep M. Trigo-Rodríguez; Joan Dergham; Maria Gritsevich; Esko Lyytinen; Elizabeth A. Silber; Iwan P. Williams
Journal: Advances in Astronomy
 (2021)

---

## Abstract

In this study, we investigate the ablation properties of bolides capable of producing meteorites. The casual dashcam recordings from many locations of the Chelyabinsk superbolide associated with the atmospheric entry of an 18 m in diameter near-Earth object (NEO) have provided an excellent opportunity to reconstruct its atmospheric trajectory, deceleration, and heliocentric orbit. In this study, we focus on the study of the ablation properties of the Chelyabinsk bolide on the basis of its deceleration and fragmentation. We explore whether meteoroids exhibiting abrupt fragmentation can be studied by analyzing segments of the trajectory that do not include a disruption episode. We apply that approach to the lower part of the trajectory of the Chelyabinsk bolide to demonstrate that the obtained parameters are consistent. To do that, we implemented a numerical (Runge–Kutta) method appropriate for deriving the ablation properties of bolides based on observations. The method was successfully tested with the cases previously published in the literature. Our model yields fits that agree with observations reasonably well. It also produces a good fit to the main observed characteristics of Chelyabinsk superbolide and provides its averaged ablation coefficientσ = 0.034 s2 km−2. Our study also explores the main implications for impact hazard, concluding that tens of meters in diameter NEOs encountering the Earth in grazing trajectories and exhibiting low geocentric velocities are penetrating deeper into the atmosphere than previously thought and, as such, are capable of producing meteorites and even damage on the ground.


---

## Body

## 1. Introduction

On February 15, 2013, our view about impact hazard was seriously challenged. While there was a sense of accomplishment for being able to forecast the close approach of 2012 DA near-Earth asteroid (NEA) within a distance of 27700 km, eventhough this NEO was discovered only one year prior, an unexpected impact with an Apollo asteroid ensued [1]. At 03:20 UTC, a superbolide, also known as the Chelyabinsk superbolide, flew over the Russian territory and Kazakhstan. The possible link between the superbolide and the 2012 DA NEA was discarded by the European Space Agency (ESA) and the Jet Propulsion Laboratory–National Aeronautics and Space Administration (JPL-NASA) from the reconstruction of the incoming fireball trajectory. The Chelyabinsk superbolide entered the atmosphere at ∼19 km/s, and according to the US sensor data (CNEOS fireball list: https://cneos.jpl.nasa.gov/fireballs/), it reached the maximum brightness at an altitude of 23.3 km with a velocity of 18.6 km/s [1, 2], also providing us with valuable samples in the form of meteorites.The existence of meteoroid streams capable of producing meteorite-dropping bolides is a hot topic in planetary science. Such streams were first proposed by Halliday [3, 4]. Their existence has important implications because they can naturally deliver to the Earth different types of rock-forming materials from potentially hazardous asteroids (PHAs). It is believed that NEOs in the Earth’s vicinity are undergoing dynamical and collisional evolution on relatively short timescales. We previously identified several NEO complexes that are producing meteorite-dropping bolides, and we hypothesize that they could have been produced during close approaches to terrestrial planets [5, 6]. Such formation scenario for this kind of asteroidal complexes is now reinforced with the recent discovery of a complex of NEOs likely associated with the NEA progenitor of the Chelyabinsk bolide [7]. The shattered pieces resulting from the disruption of NEOs visiting the inner solar system can spread along the entire parent body orbit on a time-scale of centuries [4, 8, 9]. This scenario is also consistent with the modern view of NEOs being resurfaced as consequence of close approaches [10]. Additionally, the meteorites recovered from the Chelyabinsk have a brecciated nature [11, 12] that is reminiscent of the complex collisional history and a probable rubble pile structure of the asteroid progenitor of the complex [13]. The existence of these asteroidal complexes in the near-Earth region has important implications as they could be the source of low spatial density meteoroid streams populated by large meteoroids. Such complexes could be the source of the poorly known fireball-producing radiants [14, 15]. This could have important implications for the fraction of sporadic meteoroids producing bright fireballs and in the physical mechanisms envisioned in the past [16–18].The Chelyabinsk event is also of interest because of its magnitude and energy and due to its relevance to be considered as a representative example of the most frequent outcome of the impact hazard associated with small asteroids in human timescales. Chelyabinsk also exemplifies the importance that fragmentation has for small asteroids, which can even excavate a crater on the Earth’s surface, although rarely [19–23]. Fragmentation is important as it provides a mechanism in which a significant part of the kinetic energy associated with small asteroids is released. It was certainly a very relevant process for the Tunguska event [24, 25], and in the better-known case of Chelyabinsk, most of the kinetic energy was transferred to the internal energy of the air, which is radiated as light [26].One way to study meteoroids as they enter the Earth’s atmosphere is through video observations of such events. Consequently, we are developing complementary approaches to study the dynamical behavior of video-recorded bolides in much detail. The SPanish Meteor Network (SPMN) pioneered the application of high-sensitivity cameras for detecting fireballs, and it currently maintains an online list of bright events detected over Spain, Portugal, Southern France, and Morocco since 1999 [27, 28]. For example, casual video recordings plus several still photographs of the superbolide in flight allowed us to reconstruct the heliocentric orbit of Villalbeto de la Peña meteorite in the framework of the SPMN [28]. The possibility of studying superbolides such as Chelyabinsk is a very attractive milestone to be considered. The software used in this study was developed as part of a master thesis [29] and subsequently tested and validated using several cases discussed by [29], as well as events from the 25 video and all-sky CCD stations set up over the Iberian Peninsula by the SPMN. In this context, we have been engaged in studying the dynamic behavior of meteoroids decelerating in the Earth’s atmosphere [30, 31].In this study, we study the Chelyabinsk bolide by following a Runge–Kutta method of meteor investigation similar to that developed by Bellot Rubio et al. [32]. We aim to test if that specific method is also valid for another mass range, particularly for small asteroids and large meter-sized meteoroids. We first describe our numerical model and test it against the known meteor events. We compare our code validation results to the results obtained by Bellot Rubio et al. [32] for the same dataset. We then apply our numerical model to the Chelyabinsk superbolide in order to study its dynamical behavior. For the sake of simplicity, our model considers a constant ablation coefficient and the shape factor, even though these parameters could vary in different ablation stages [33–35].This study is structured as follows: the data reduction and the theoretical approach pertaining to the Chelyabinsk bolide are described in the next section. InSection 3, the main implications of this work in the context of fireballs, meteorites, and NEO research are discussed. We use the model to determine the fireball flight parameters, and by studying the deceleration, we also obtain the ablation coefficient. Finally, the conclusions of this work are presented in Section 4.


## 2. Data Reduction, Theoretical Approach, and Observations

The Chelyabinsk superbolide was an unexpected daylight superbolide as many other unpredicted meteorite-dropping bolides in history. Fortunately, numerous casual video recordings of the bolide trajectory from the ground were obtained, given the nowadays common dashcams available in private motor vehicles in Russia. According to the video recordings available, it is possible to study the atmospheric trajectory and deceleration carefully, allowing the reconstruction of the heliocentric orbit in record time [2, 26].


### 2.1. Single Body Theory

There are two main approaches in the study of the dynamic properties of meteors during atmospheric interaction, the quasicontinuous fragmentation (QCF) theory introduced by Novikov et al. [37], which was later extended by Babadzhanov [38], and the single body theory described by Bronshsten [39]. There have been discrepancies in the applicability conditions of both methods: the single body works with basic differential equations, while the QCF uses semiempirical formulas studying only the luminosity produced by the meteor. The main difference is that the single body theory obtains smaller dynamical masses than the QCF method. Thus far, neither approach is prevalent, and the reason why the theories do not converge could be attributed to the contribution of other key processes such as the fragmentation and deceleration of meteoroids during ablation or to the poorly constrained values of the bulk density and/or the luminous efficiency coefficient [40–42].We remark that the initial dynamic mass estimate or the preatmospheric size can be derived using the methods described in other works [43–48]; hence, we leave it out of the scope of the present model. We also note that alternative models accounting for ablation have been recently developed; however, further discussion on that topic is beyond the scope of this study, and the reader is referred to the following literature [46, 47, 49–52]. As noted in the introduction, high-strength meteoroids from asteroids or planetary bodies exhibit a very different behavior than fragile dust aggregates originating from comets [53–57].


### 2.2. The Role of Fragmentation

The fragmentation of meteoroids was studied in detail by various authors [38, 40]. After analyzing different photographic observations, Levin [40] distinguished four possible types of fragmentation: (a) the decay of the meteoroid into large nonfragmenting debris, (b) the progressive disintegration of the original meteoroid into fragments that continue to crumble into smaller fragments, (c) the instantaneous ejection of a large number of small particles, which, when affecting the entire meteoroid, is called a catastrophic disruption, and finally, (d) the quasicontinuous fragmentation which consists of the gradual release of a large number of small particles from the surface and their subsequent evaporation due to high temperatures associated with the shock thermal wave formed around the body.In practice, a combination of two or more types of fragmentation can be observed in a given meteor event. In fact, it is possible to observe that (a) and (c) fragmentation types described in the preceding paragraph could occur more than once for the same meteor event. The analysis of meteors performed by Jacchia [58] using Super-Schmidt cameras demonstrated that the single body theory does not work for cases which suffer abrupt types of fragmentation along the trajectory. As a direct consequence, meteoroids exhibiting fragmentation of the first (a), second (b), and third (c) kind should not be studied using this simplified single body theory. When a meteoroid suffers an abrupt fragmentation, the main body loses mass instantaneously, and therefore, the single body equations cannot be applied since the condition of continuity of mass is not satisfied. Thus, the cases with possible abrupt fragmentation episode(s) will not be examined in this work.


### 2.3. The Single Body Theory: The Drag and Mass Loss Equations

The dynamic behavior of a meteoroid as it interacts with the Earth’s atmosphere is described using the drag and mass loss equations. These equations, as presented by Bronshten [39], are as follows:(1)dvdt=−Kρairm−1/3·v2,(2)dmdt=−σKρairm2/3·v3.where K is the shape-density coefficient, ρair is the air density, m is the meteoroid mass, v is the instantaneous velocity, and σ is the ablation coefficient.By using equations (1) and (2), the parameters to be identified are K and σ. The ablation coefficient defines the loss of mass for the bolide as it penetrates the atmosphere; the larger the value is, the more the mass will be ablated for a given velocity. The value of the ablation coefficient depends on various factors and is expressed as(3)σ=Λ2ΓQ,where Λ is the heat transfer coefficient, Γ is the drag coefficient, and Q is the heat of ablation.The shape-density coefficient depends on the shape and density of the meteoroid and is expressed as(4)K=sAρ2/3,where A is the shape factor, s is the cross-section area, and ρ is the meteoroid bulk density.We must point out that the observational data obtained from the reconstruction of the trajectories of meteors using CCD or video cameras is basically the frame to frame speed of the bolide as a function of the height, requiring another equation to link the time with the altitude:(5)dhdt=−v·cosz,where z is the zenith angle.By substituting equation (5) into equations (1) and (2), the following expressions are obtained:(6)dvdh=Kρairm−1/3·vcosz,(7)dmdh=σKρairm2/3·v2cosz.Next, by dividing equation (2) by equation (1), we obtain(8)dmdv=mσv.Solving this differential equation with the boundary condition ofm=mo, when v=vo, the following is obtained:(9)m=mo·e−1/2σvo2−v2.Now, we combine equations (9) and (6) to derive(10)dvdh=Kρairmo·e−1/2σvo2−v2−1/3⋅v⋅1cosz.In order to obtain the value ofK and σ, we use equation (10), as it uses substitution for the dependence on the instantaneous mass, and thus, we deal with one equation instead of two. Furthermore, there is the initial mass, which is an important parameter to be studied. Equation (10) directly links the deceleration of the meteoroid as a function of different parameters, in particular the zenith angle. The last term is of particular interest as it modulates the full equation. For the vertical entry (z = 0°), the deceleration is maximized, while for z close to 90°, the deceleration is minimized. For that reason, large meteoroids at grazing angles are able to follow extremely long trajectories or even escape again into space, such as the Grand Tetons superbolide that spent nearly two minutes traveling over several states of the USA and Canada on August 10, 1972 [11, 59].However, by using the concepts introduced above, it is not possible to obtain the values of initial mass (mo) and K. The parameters that can be found are only the expressions of mo−1/3·K and σ. Consequently, another equation is required to obtain K and mo separately. The remaining expression is the photometric equation:(11)I=τ2−dmdt·v2,where τ is the luminous efficiency. This equation is assuming v = constant, and it is often applied to small meteors. The luminous efficiency is obtained empirically; thus, any deviations in that value might produce significant changes in the results. This equation will be referred to later as a part of the mass determination. Now, we focus only on equation (10).We define theK′ as(12)K′=K⋅mo−1/3.Then, the equation to work with is(13)dvdh=ρairK′e−1/2σvo2−v2−1/3⋅v⋅1cosz,with variables K′ and σ.


## 2.1. Single Body Theory

There are two main approaches in the study of the dynamic properties of meteors during atmospheric interaction, the quasicontinuous fragmentation (QCF) theory introduced by Novikov et al. [37], which was later extended by Babadzhanov [38], and the single body theory described by Bronshsten [39]. There have been discrepancies in the applicability conditions of both methods: the single body works with basic differential equations, while the QCF uses semiempirical formulas studying only the luminosity produced by the meteor. The main difference is that the single body theory obtains smaller dynamical masses than the QCF method. Thus far, neither approach is prevalent, and the reason why the theories do not converge could be attributed to the contribution of other key processes such as the fragmentation and deceleration of meteoroids during ablation or to the poorly constrained values of the bulk density and/or the luminous efficiency coefficient [40–42].We remark that the initial dynamic mass estimate or the preatmospheric size can be derived using the methods described in other works [43–48]; hence, we leave it out of the scope of the present model. We also note that alternative models accounting for ablation have been recently developed; however, further discussion on that topic is beyond the scope of this study, and the reader is referred to the following literature [46, 47, 49–52]. As noted in the introduction, high-strength meteoroids from asteroids or planetary bodies exhibit a very different behavior than fragile dust aggregates originating from comets [53–57].


## 2.2. The Role of Fragmentation

The fragmentation of meteoroids was studied in detail by various authors [38, 40]. After analyzing different photographic observations, Levin [40] distinguished four possible types of fragmentation: (a) the decay of the meteoroid into large nonfragmenting debris, (b) the progressive disintegration of the original meteoroid into fragments that continue to crumble into smaller fragments, (c) the instantaneous ejection of a large number of small particles, which, when affecting the entire meteoroid, is called a catastrophic disruption, and finally, (d) the quasicontinuous fragmentation which consists of the gradual release of a large number of small particles from the surface and their subsequent evaporation due to high temperatures associated with the shock thermal wave formed around the body.In practice, a combination of two or more types of fragmentation can be observed in a given meteor event. In fact, it is possible to observe that (a) and (c) fragmentation types described in the preceding paragraph could occur more than once for the same meteor event. The analysis of meteors performed by Jacchia [58] using Super-Schmidt cameras demonstrated that the single body theory does not work for cases which suffer abrupt types of fragmentation along the trajectory. As a direct consequence, meteoroids exhibiting fragmentation of the first (a), second (b), and third (c) kind should not be studied using this simplified single body theory. When a meteoroid suffers an abrupt fragmentation, the main body loses mass instantaneously, and therefore, the single body equations cannot be applied since the condition of continuity of mass is not satisfied. Thus, the cases with possible abrupt fragmentation episode(s) will not be examined in this work.


## 2.3. The Single Body Theory: The Drag and Mass Loss Equations

The dynamic behavior of a meteoroid as it interacts with the Earth’s atmosphere is described using the drag and mass loss equations. These equations, as presented by Bronshten [39], are as follows:(1)dvdt=−Kρairm−1/3·v2,(2)dmdt=−σKρairm2/3·v3.where K is the shape-density coefficient, ρair is the air density, m is the meteoroid mass, v is the instantaneous velocity, and σ is the ablation coefficient.By using equations (1) and (2), the parameters to be identified are K and σ. The ablation coefficient defines the loss of mass for the bolide as it penetrates the atmosphere; the larger the value is, the more the mass will be ablated for a given velocity. The value of the ablation coefficient depends on various factors and is expressed as(3)σ=Λ2ΓQ,where Λ is the heat transfer coefficient, Γ is the drag coefficient, and Q is the heat of ablation.The shape-density coefficient depends on the shape and density of the meteoroid and is expressed as(4)K=sAρ2/3,where A is the shape factor, s is the cross-section area, and ρ is the meteoroid bulk density.We must point out that the observational data obtained from the reconstruction of the trajectories of meteors using CCD or video cameras is basically the frame to frame speed of the bolide as a function of the height, requiring another equation to link the time with the altitude:(5)dhdt=−v·cosz,where z is the zenith angle.By substituting equation (5) into equations (1) and (2), the following expressions are obtained:(6)dvdh=Kρairm−1/3·vcosz,(7)dmdh=σKρairm2/3·v2cosz.Next, by dividing equation (2) by equation (1), we obtain(8)dmdv=mσv.Solving this differential equation with the boundary condition ofm=mo, when v=vo, the following is obtained:(9)m=mo·e−1/2σvo2−v2.Now, we combine equations (9) and (6) to derive(10)dvdh=Kρairmo·e−1/2σvo2−v2−1/3⋅v⋅1cosz.In order to obtain the value ofK and σ, we use equation (10), as it uses substitution for the dependence on the instantaneous mass, and thus, we deal with one equation instead of two. Furthermore, there is the initial mass, which is an important parameter to be studied. Equation (10) directly links the deceleration of the meteoroid as a function of different parameters, in particular the zenith angle. The last term is of particular interest as it modulates the full equation. For the vertical entry (z = 0°), the deceleration is maximized, while for z close to 90°, the deceleration is minimized. For that reason, large meteoroids at grazing angles are able to follow extremely long trajectories or even escape again into space, such as the Grand Tetons superbolide that spent nearly two minutes traveling over several states of the USA and Canada on August 10, 1972 [11, 59].However, by using the concepts introduced above, it is not possible to obtain the values of initial mass (mo) and K. The parameters that can be found are only the expressions of mo−1/3·K and σ. Consequently, another equation is required to obtain K and mo separately. The remaining expression is the photometric equation:(11)I=τ2−dmdt·v2,where τ is the luminous efficiency. This equation is assuming v = constant, and it is often applied to small meteors. The luminous efficiency is obtained empirically; thus, any deviations in that value might produce significant changes in the results. This equation will be referred to later as a part of the mass determination. Now, we focus only on equation (10).We define theK′ as(12)K′=K⋅mo−1/3.Then, the equation to work with is(13)dvdh=ρairK′e−1/2σvo2−v2−1/3⋅v⋅1cosz,with variables K′ and σ.


## 3. Results and Discussion

### 3.1. Numerical Approximation

In this section, we develop a numerical approximation with the aim to describe the meteoroid flight in the atmosphere. Our goal is to obtain the solution that can be used to better understand that physical process. Subsequently, our intention is to develop a numerical approach that can be very valuable in predicting the variation of the parameters along the trajectory segments versus analytical “whole-trajectory smoothing.”Equation (13) is the expression used to find the physical parameters. To test our model, we use meteor data pertaining to the meteoroid velocity at different altitudes taken from the literature. In this regard, the accurate trajectory/velocity meteor data obtained via high-resolution Super-Schmidt cameras [58] can be utilized for a case study. These meteor trajectory data are used to optimize the procedure to obtain the values of K′ and σ that fit the best the actual data points. The procedure presented here produces many synthetic curves. Subsequently, in order to optimize the computing time, a solution to make the equation converge to the real data is found.


#### 3.1.1. The Runge–Kutta Method Implementation

The Runge–Kutta method is an iterative technique for the approximation and solution of ordinary differential equations. The method was first developed by Runge [60] and Kutta [61].The Runge–Kutta approximation provides a solution at a determined point of altitude. Application of the Runge–Kutta method requires the initial conditions to be known:(14)y′=ft,y⋅;yo;to.In our case, the initial conditions will be the initial velocity and the altitude of the bolide when ablation starts, “synthetically” written as(15)dvdh=fh,v;vho;ho˙.We then need to choose a step size (p) that should be comparable with the data resolution to allow a better comparison between the model and the observations. The step size defines how many integration steps need to be implemented before the final solution is reached. The smaller the step size, the more accurate the solution will be, noting that this will also increase the computing time. The step size corresponds to the characteristic Δh that according to Bellot Rubio et al. [32] can be chosen around a few hundred of meters, 100–300 m.Once the step size is defined, we define the model coefficients as follows:(16)l1=ftn,yn,l2=ftn+12p,yn+p12l1,l3=ftn+12p,yn+p12l2,l4=ftn+p,yn+pl3.For our case, the function to be studied is(17)fhn,vn=K′ρhe−1/2σvo2−v2−1/3⋅v⋅1cosz,and the coefficients are calculated as(18)l1=fhn,vn,l2=fhn+12p,vn+p12l1,l3=fhn+12p,vn+p12l2,l4=fhn+p,vn+pl3.Once the coefficients (equation (16)) are calculated, we compute the solution for the point yn+1 using the following formula:(19)yn+1=yn+16pl1+2l2+2l3+l4.For our case (equation (18)), this is(20)vn+1=vn+16pl1+2l2+2l3+l4.The result obtained is the solution for the point (hn+1, vn+1), which becomes the initial condition to find the numerical approximation for the next point. The procedure is repeated until the desired value is reached.


#### 3.1.2. Validating the Code Using the Jacchia Catalog

Once the procedure is defined, we need to validate the code by comparing the results with the previously published data. We use a catalog of very precise photographic trajectories of meteors [58], hereafter referred to as the JVB catalog, obtained using high spatial resolution Super-Schmidt cameras. Jacchia’ [58] study synthesizes the physical inferred parameters for 413 meteors between −5 and +2, 5 stellar magnitude ranges obtained during a multistation meteor network operated in the fifties and sixties in New Mexico, USA. The data provide the meteor velocity and magnitude as a function of the altitude, the derived preatmospheric velocity, the deceleration, and some additional information for the observed meteors. All JVB catalog events were named using a number. In this project, we have used the same Jacchia and Whipple [58] numbering, but included a J at the beginning for the computational reasons. For example, we later discuss meteor J8945 (listed in JVB as 8945).Equation (13) also requires the knowledge of air density. We adopted a general model widely used in meteor studies [45], the United States (US) standard atmosphere [62]. The US standard atmosphere was originally developed in 1958 by the US Committee on Extension to the Standard Atmosphere and improved in 1976. It is a series of tables which approximate the values for atmospheric temperature, density, pressure, and other properties over a wide range of altitudes.


#### 3.1.3. Autofitting Procedure

We have defined a way to transform the differential equation into an expression that can be iteratively computed. The aim is to find a result for whichK′ and σ produce the closest curve fitting to the observational data points. In order to find the best matching values, we introduce the following autofitting procedure. We start by choosing two random values of K′ and σ set as initial approximation. The meteor data to be investigated include the velocity of a meteor at specified points of altitude; we compare this velocity with the velocity simulated by our code at the same altitude points. The error factor is then calculated as follows:(21)ξ=vd−vfK,σ2,where vd is the velocity inferred from the measured data, and vf is the computed velocity.Furthermore, we introduce the increment factors forK′ and σ, which are defined as ΔK and Δσ. Following the same procedure as earlier, we compute the error factor for(22)ξ=vd−vfK+ΔK,σ2ξ=vd−vfK−ΔK,σ2,ξ=vd−vfK,σ+Δσ2ξ=vd−vfK,σ−Δσ2,ξ=vd−vfK+ΔK,σ+Δσ2ξ=vd−vfK−ΔK,σ−Δσ2,ξ=vd−vfK+ΔK,σ−Δσ2ξ=vd−vfK−ΔK,σ+Δσ2.In principle, we created a 2D matrix of errors. The error parameter can be shown in a table for better visualization of the algorithm (Figure1). Once all the error values are computed, the search for the minimum value is performed, and the result is considered the new centered value for the next computation iteration.Figure 1
Basic scheme of the error matrix.We repeat the same procedure for the next centered value, until we reach the point where the minimum will be centered in the middle of the matrix. Consequently, the minimum error value will correspond to the sought values ofK′ and σ. If very small increments of K′ and σ are used, the solution will be more accurate (at the cost of increased computing time). If we set large increments of K′ and σ, the solution will be reached faster, but at the expense of the resolution of the solution. In order to deal with this conundrum, the code is optimized such that it is set to work with large increments of K′ and σ at the beginning. Once the ‘first approximation’ solution is found, the code switches to smaller increments until the optimal resolution is reached.Figure2 shows an example of the autofitting procedure. We have chosen the J8945 meteor for the comparison with published data [32]. Other cases were studied and compared in Dergham [29]. By comparing the velocity vs. altitude graphs, it is evident that the fitting curves are very similar to the observed data.Figure 2
Normalized velocity vs. altitude fit for the meteor J8945. The result is a good fit to the previously published Figure1 of Bellot Rubio et al. [32].


#### 3.1.4. Abrupt Disruption Cases

We have presented a model capable of obtaining some parameters for meteoroids. However, as previously mentioned, not all meteoroids can be studied using this particular model because if they undergo fragmentation, the results might be skewed. Bellot Rubio et al. [32] also mentioned that there are a significant number of cases in the JVB catalog that cannot be fitted, likely as a consequence of undergoing abrupt disruptions. Figure 3 shows the results of the velocity curve using our model for the case J4141. The results are also compared with that obtained by Bellot Rubio et al. [32].Figure 3
Comparison of results for the meteoroid J4141. The model cannot produce a well-fitting solution since this meteoroid likely suffered an abrupt disruption. This figure can be compared with that shown in Figure1 of Bellot Rubio et al. [32].Despite the difficulties to obtain well-fitting solutions for some events, it is remarkable that our model is able to identify and produce solutions for the events undergoing quasicontinuous fragmentation. A possible way to study meteoroids with abrupt fragmentation is to focus on different segments of the trajectory that do not include a disruption episode. We will apply that approach to the lower part of the trajectory of the Chelyabinsk bolide to demonstrate that the obtained parameters are consistent.


### 3.2. The Ablation Coefficient

Equation (13) has several unknowns. The ability of a meteoroid to produce light can be linked to its mass loss or the ablation coefficient σ [34, 42]. In principle, the ablation coefficient reflects how rapidly the meteoroid loses its mass as it interacts with the atmosphere. Low values of the ablation coefficient indicate that the object is losing a lesser amount of mass compared to that having a higher value of the ablation coefficient. The ablation coefficient is normally represented in units of s2 km−2 and can also be expressed through the dimensionless mass loss parameter [44, 63] used through the range of meteor physics studies. The value of the ablation coefficient depends of many factors, such as chemical composition, grain size, density, porosity, and body shape, among others. In general, the values of the ablation coefficient range between 0.01 and 0.3 s2 km−2 [25]. In order to exemplify this, we applied different ablation coefficients to a 1 g meteoroid with the preatmospheric velocity of 25 km/s and which starts to decelerate at the altitude of 100 km. The results are shown in Figure 4.Figure 4
Comparison of meteoroids having the same initial conditions (m0 = 1 g, v0 = 25 km/s) but different ablation coefficients.In general, the larger the ablation coefficient is, the faster the body decelerates due to a more rapid mass loss. Consequently, the mass of the body decreases due to ablation; this is described by using the ablation coefficient, and the drag force imposed by the atmosphere has a larger effect. Table1 shows the comparison of our results to several events described in the JVB catalog.Table 1
Comparison of the ablation coefficients for various meteor cases from the JVB catalog.
Meteor IDJacchia (s2·km−2)This study (s2·km−2)J68820.08120.075J69590.03310.0382J72160.05010.079J89450.03540.036J90300.05490.0542J71610.03540.0381


### 3.3. Deceleration, Normalized Instantaneous Mass, and Mass Loss Rate

In this section, the effect of deceleration is examined in more detail. Given that we have the velocity as a function of altitude along the trajectory, we can study deceleration from the behavior of the velocity curve. This is particularly useful as many meteor processing algorithms and detection methods provide velocity values sequentially [15, 64]. At a certain point hi, the code computes an increment of velocity over the increment of distance at the points immediately before and after. This can be expressed as(23)dvdhhi≅ΔvΔh=vh2−vh1h2−h1.Considering all the points with known velocity and altitude along the trajectory as input data, equation (23) can be applied directly. Figure 5(a) shows the velocity curve as a function of altitude for the J8945 meteoroid.Figure 5
Curves of the acceleration (a), mass evolution (b), and relative mass loss rate (c) as a function of altitude for the meteoroid J8945.
(a)(b)(c)The normalized instantaneous mass (m/m0) is the next quantity to be studied. The expression for the normalized instantaneous mass can be derived by rearranging equation (9):(24)mmo=e−1/2σvo2−v2.Equation (24) expresses the normalized instantaneous mass as a function of velocity, whereas the values of velocity are the functions of altitude. Figure 5(b) shows the behavior of normalized mass expressed as a function of altitude for the J8945 event.We define the normalized mass loss rate as the derivative of the relative mass over the altitude. This value is computed as follows:(25)dm/moh=Δm/moΔh=m/moh2−m/moh1h2−h1.Figure5(c) shows the direct application of this approach for the J8945 event.


### 3.4. Application to the Chelyabinsk Superbolide

We apply the Runge–Kutta code developed in this work to the famous Chelyabinsk superbolide. On February 15, 2013, it was predicted that the NEA-2012 DA14, discovered a year prior by the Observatory Astronòmic de Mallorca, would approach the Earth within a close minimum distance of only 27700 km. However, while all the attention was focused on anticipating that encounter, another NEA unexpectedly entered the atmosphere over central Asia on February 15, 2013 at 03:20:33 UTC. The bolide disintegrated in the proximity of the city of Chelyabinsk [1]. The Chelyabinsk bolide reached the brightness magnitude of −28, which is brighter than the moon (Figure 6).Figure 6
A casual photograph of the Chelyabinsk bolide taken by Marat Ahmetvaleev. (a) The abrupt fragmentation at a height of 23 km is clearly distinguishable by the main flare. (b) Pictures of the dust trail by the same photographer reveal multiple paths produced by the fragments. Both pictures are courtesy of the author.
(a)(b)As the days passed and the orbits were calculated, scientists discarded a possible association between the two NEAs, as they presented very different heliocentric orbits. Thanks to video cameras (dash cams) found in the majority of Russian motor vehicles and surveillance cameras placed on buildings, the initial trajectory of the bolide was reconstructed and the orbit was determined [2].After the superbolide sightings, many people uploaded various videos to Internet. Since the geographical location of the recorded videos was known, we reconstructed the bolide trajectory, obtaining the values of velocity as a function of altitude. As shown in Table2, the data were obtained from the analysis of our video compilation, the bolide velocity along the terminal part of its trajectory just after the massive fragmentation event taking place at the altitude of 26 km. The step size was determined by the video frame rate, corresponding to differences in height of 200–150 m. Table 2 shows these data.Table 2
Dynamic data of the Chelyabinsk superbolide.
Height (km)V (km/s)Height (km)v (km/s)18.9814.0415.669.7318.7813.8615.539.4618.5813.6815.399.2018.3813.4915.268.9418.1813.2915.138.6817.9913.0915.018.4217.8012.8814.898.1717.6212.6614.777.9217.4412.4414.667.6717.2612.2214.557.4317.0811.9914.457.1916.9111.7514.346.9616.7411.5114.246.7416.5811.2714.156.5216.4211.0214.066.3116.2610.7613.976.1116.1010.5113.885.9215.9510.2513.805.73Figure7(a) depicts the dynamical data after the disruption had occurred at an altitude of 26 km, and the fit is obtained by our model. The graph shows a quite uniform behavior of the Chelyabinsk superbolide after the main fragmentation event.Figure 7
Velocity evolution (a), mass evolution (b), and mass loss rate (c) as a function of altitude for the Chelyabinsk superbolide.
(a)(b)(c)By studying the dynamic curve, the ablation coefficient can be obtained. The derived value is(26)σ=0.034s2⋅km−2.It is quite remarkable that this value, derived for the lower trajectory, provides a similar ablation coefficient as that for the fireballs at much higher altitudes, even though Chelyabinsk was the deepest penetrating bolide ever documented, still emitting light even as it reached the troposphere. Notably, the atmospheric density in these lower regions is about four orders of magnitude higher.The normalized mass evolution for that lower part of the trajectory is shown in Figure7(b), and the mass loss rate evolution is shown in Figure 7(c). In particular, it is certainly encouraging that the model predicts quite well the ablation behavior of the Chelyabinsk bolide in the lower part of its atmospheric trajectory. The light curve of the Chelyabinsk event was normalized using the US government sensor data at a peak of brightness value of 2.7 we1013 W s·r−1, corresponding to an absolute astronomical magnitude of −28 [1]. According to the NASA JPL Chelyabinsk final report [65], the maximum brightness was achieved at an altitude of 23.3 km. This is consistent with our results as the maximum inferred value of mass loss rate from our dynamic data that occur at an altitude of ∼23.5 km (Figure 7(c)).It is well known that the maximum brightness is achieved shortly after the meteoroid catastrophically breaks apart due to the fragmented/pulverized material being exposed to the heat generated by the resulting shock wave.An interesting conclusion that can be obtained directly from these results is the importance of the atmosphere. As mentioned before, the faster the meteoroid, the more rapid the ablation process. Thus, the atmosphere could effectively shield the Earth from very fast impacts since such objects are preferentially and more quickly ablated. However, the less favorable cases are very large objects (especially if their preatmospheric velocity is low), such as the Tunguska impactor that produced an airburst over Siberia in 1908 when it entered the atmosphere at a velocity of ∼30 km/s [66]. The Chelyabinsk superbolide had the preatmospheric velocity of 19 km/s, and its trajectory was consistent with a grazing geometry (high zenith angle).To exemplify this, we have plotted the entry of the Chelyabinsk superbolide for different initial velocities (Figure8). Of course, disrupting NEOs can produce fragments tens of meters in size. If these fragments encounter our planet under certain geometric conditions, they might become a significant source of damage on the ground and casualties. Thus, identifying the existence of asteroidal complexes in the near-Earth environment is crucial for better assessment of impact hazard.Figure 8
Simulation of the Chelyabinsk superbolide with different initial velocities (v = 17.1, 23, and 27 km/s).If the Chelyabinsk superbolide entered at a higher velocity, it would have slowed down faster due to a more rapid mass loss. According to our model, the maximum brightness of a meteoroid occurs when the mass loss reaches the peak. Consequently, the model predicts that a Chelyabinsk-like asteroid moving at a lower velocity could have more destructive potential on the Earth’s surface (Figure8). On the other hand, a projectile of the same rocky composition and moving at higher velocity could have undergone a much faster ablation process, culminating in an explosive airburst, similar to that occurring over Tunguska. Such conclusions imply that the efficiency of the Earth’s atmosphere to shield us from dangerous tens of meters in size asteroids that strongly depend on the relative velocity of the encounter with our planet.It is important to remark that new improvements in detection of fireballs from space could provide an additional progress in studying the luminous efficiency of bolides [34]. In fact, several detections of the Chelyabinsk bolide from space have been reported [67, 68]. Future studies of common events detected from both the ground and space could constrain the role of the observing geometry in signal loss and possible biases in the subsequent determination of velocity, radiated energy, and determination of the orbital elements.


### 3.5. Implications for Impact Hazard

Considering the destructive nature of large extraterrestrial objects (e.g., [24]), it is crucial to identify the existence of asteroidal complexes in the near-Earth environment. Indeed, the Chelyabinsk scale or larger extraterrestrial bodies (10s of meters in size) are of special interest to the planetary defense community because depending on their orbital geometry and impact angle, these objects could pose significant hazard for humans and infrastructure on the ground. Thus, studying the near-Earth environment along with the well-documented events such as the Chelyabinsk superbolide can shed more light on the dynamic processes, as well as fundamental properties of these objects.From the study of meteorite falls and the relative absence of small impact craters [19], we know that meter-sized asteroids are efficiently fragmented when they penetrate into the stratosphere at hypervelocity. The loading pressure in front of the body produces the rock fracture when it overcomes its tensile strength. As a natural consequence, meteorite falls often deliver tens to hundreds of stones, just after this type of break-ups [69, 70]. We have previously described such behavior when discussing the evolution of Chelyabinsk. Additionally, there is a consensual view that the terrestrial atmosphere behaves as an efficient shield for meter to 10s of meters in diameter projectiles. Despite of this, we still need to increase our statistics as more rarely, and probably for low-velocity projectiles under favorable geometric circumstances, crater excavation is still possible. A good example of this was the so-called Carancas event: an impact crater excavated in the Perú’s Altiplane by a one meter-sized chondritic meteorite [71]. Being a quite unusual impact event, it is also likely that a significant number of these events are rarely studied because they happen in remote locations and remain unnoticed. It is obvious that Chelyabinsk and other well-recorded events (e.g., by fireball networks) provide an opportunity to understand in the required detail the behavior of meter-sized meteoroids and their capacity to directly cause human injuries and even casualties.


## 3.1. Numerical Approximation

In this section, we develop a numerical approximation with the aim to describe the meteoroid flight in the atmosphere. Our goal is to obtain the solution that can be used to better understand that physical process. Subsequently, our intention is to develop a numerical approach that can be very valuable in predicting the variation of the parameters along the trajectory segments versus analytical “whole-trajectory smoothing.”Equation (13) is the expression used to find the physical parameters. To test our model, we use meteor data pertaining to the meteoroid velocity at different altitudes taken from the literature. In this regard, the accurate trajectory/velocity meteor data obtained via high-resolution Super-Schmidt cameras [58] can be utilized for a case study. These meteor trajectory data are used to optimize the procedure to obtain the values of K′ and σ that fit the best the actual data points. The procedure presented here produces many synthetic curves. Subsequently, in order to optimize the computing time, a solution to make the equation converge to the real data is found.


### 3.1.1. The Runge–Kutta Method Implementation

The Runge–Kutta method is an iterative technique for the approximation and solution of ordinary differential equations. The method was first developed by Runge [60] and Kutta [61].The Runge–Kutta approximation provides a solution at a determined point of altitude. Application of the Runge–Kutta method requires the initial conditions to be known:(14)y′=ft,y⋅;yo;to.In our case, the initial conditions will be the initial velocity and the altitude of the bolide when ablation starts, “synthetically” written as(15)dvdh=fh,v;vho;ho˙.We then need to choose a step size (p) that should be comparable with the data resolution to allow a better comparison between the model and the observations. The step size defines how many integration steps need to be implemented before the final solution is reached. The smaller the step size, the more accurate the solution will be, noting that this will also increase the computing time. The step size corresponds to the characteristic Δh that according to Bellot Rubio et al. [32] can be chosen around a few hundred of meters, 100–300 m.Once the step size is defined, we define the model coefficients as follows:(16)l1=ftn,yn,l2=ftn+12p,yn+p12l1,l3=ftn+12p,yn+p12l2,l4=ftn+p,yn+pl3.For our case, the function to be studied is(17)fhn,vn=K′ρhe−1/2σvo2−v2−1/3⋅v⋅1cosz,and the coefficients are calculated as(18)l1=fhn,vn,l2=fhn+12p,vn+p12l1,l3=fhn+12p,vn+p12l2,l4=fhn+p,vn+pl3.Once the coefficients (equation (16)) are calculated, we compute the solution for the point yn+1 using the following formula:(19)yn+1=yn+16pl1+2l2+2l3+l4.For our case (equation (18)), this is(20)vn+1=vn+16pl1+2l2+2l3+l4.The result obtained is the solution for the point (hn+1, vn+1), which becomes the initial condition to find the numerical approximation for the next point. The procedure is repeated until the desired value is reached.


### 3.1.2. Validating the Code Using the Jacchia Catalog

Once the procedure is defined, we need to validate the code by comparing the results with the previously published data. We use a catalog of very precise photographic trajectories of meteors [58], hereafter referred to as the JVB catalog, obtained using high spatial resolution Super-Schmidt cameras. Jacchia’ [58] study synthesizes the physical inferred parameters for 413 meteors between −5 and +2, 5 stellar magnitude ranges obtained during a multistation meteor network operated in the fifties and sixties in New Mexico, USA. The data provide the meteor velocity and magnitude as a function of the altitude, the derived preatmospheric velocity, the deceleration, and some additional information for the observed meteors. All JVB catalog events were named using a number. In this project, we have used the same Jacchia and Whipple [58] numbering, but included a J at the beginning for the computational reasons. For example, we later discuss meteor J8945 (listed in JVB as 8945).Equation (13) also requires the knowledge of air density. We adopted a general model widely used in meteor studies [45], the United States (US) standard atmosphere [62]. The US standard atmosphere was originally developed in 1958 by the US Committee on Extension to the Standard Atmosphere and improved in 1976. It is a series of tables which approximate the values for atmospheric temperature, density, pressure, and other properties over a wide range of altitudes.


### 3.1.3. Autofitting Procedure

We have defined a way to transform the differential equation into an expression that can be iteratively computed. The aim is to find a result for whichK′ and σ produce the closest curve fitting to the observational data points. In order to find the best matching values, we introduce the following autofitting procedure. We start by choosing two random values of K′ and σ set as initial approximation. The meteor data to be investigated include the velocity of a meteor at specified points of altitude; we compare this velocity with the velocity simulated by our code at the same altitude points. The error factor is then calculated as follows:(21)ξ=vd−vfK,σ2,where vd is the velocity inferred from the measured data, and vf is the computed velocity.Furthermore, we introduce the increment factors forK′ and σ, which are defined as ΔK and Δσ. Following the same procedure as earlier, we compute the error factor for(22)ξ=vd−vfK+ΔK,σ2ξ=vd−vfK−ΔK,σ2,ξ=vd−vfK,σ+Δσ2ξ=vd−vfK,σ−Δσ2,ξ=vd−vfK+ΔK,σ+Δσ2ξ=vd−vfK−ΔK,σ−Δσ2,ξ=vd−vfK+ΔK,σ−Δσ2ξ=vd−vfK−ΔK,σ+Δσ2.In principle, we created a 2D matrix of errors. The error parameter can be shown in a table for better visualization of the algorithm (Figure1). Once all the error values are computed, the search for the minimum value is performed, and the result is considered the new centered value for the next computation iteration.Figure 1
Basic scheme of the error matrix.We repeat the same procedure for the next centered value, until we reach the point where the minimum will be centered in the middle of the matrix. Consequently, the minimum error value will correspond to the sought values ofK′ and σ. If very small increments of K′ and σ are used, the solution will be more accurate (at the cost of increased computing time). If we set large increments of K′ and σ, the solution will be reached faster, but at the expense of the resolution of the solution. In order to deal with this conundrum, the code is optimized such that it is set to work with large increments of K′ and σ at the beginning. Once the ‘first approximation’ solution is found, the code switches to smaller increments until the optimal resolution is reached.Figure2 shows an example of the autofitting procedure. We have chosen the J8945 meteor for the comparison with published data [32]. Other cases were studied and compared in Dergham [29]. By comparing the velocity vs. altitude graphs, it is evident that the fitting curves are very similar to the observed data.Figure 2
Normalized velocity vs. altitude fit for the meteor J8945. The result is a good fit to the previously published Figure1 of Bellot Rubio et al. [32].


### 3.1.4. Abrupt Disruption Cases

We have presented a model capable of obtaining some parameters for meteoroids. However, as previously mentioned, not all meteoroids can be studied using this particular model because if they undergo fragmentation, the results might be skewed. Bellot Rubio et al. [32] also mentioned that there are a significant number of cases in the JVB catalog that cannot be fitted, likely as a consequence of undergoing abrupt disruptions. Figure 3 shows the results of the velocity curve using our model for the case J4141. The results are also compared with that obtained by Bellot Rubio et al. [32].Figure 3
Comparison of results for the meteoroid J4141. The model cannot produce a well-fitting solution since this meteoroid likely suffered an abrupt disruption. This figure can be compared with that shown in Figure1 of Bellot Rubio et al. [32].Despite the difficulties to obtain well-fitting solutions for some events, it is remarkable that our model is able to identify and produce solutions for the events undergoing quasicontinuous fragmentation. A possible way to study meteoroids with abrupt fragmentation is to focus on different segments of the trajectory that do not include a disruption episode. We will apply that approach to the lower part of the trajectory of the Chelyabinsk bolide to demonstrate that the obtained parameters are consistent.


## 3.1.1. The Runge–Kutta Method Implementation

The Runge–Kutta method is an iterative technique for the approximation and solution of ordinary differential equations. The method was first developed by Runge [60] and Kutta [61].The Runge–Kutta approximation provides a solution at a determined point of altitude. Application of the Runge–Kutta method requires the initial conditions to be known:(14)y′=ft,y⋅;yo;to.In our case, the initial conditions will be the initial velocity and the altitude of the bolide when ablation starts, “synthetically” written as(15)dvdh=fh,v;vho;ho˙.We then need to choose a step size (p) that should be comparable with the data resolution to allow a better comparison between the model and the observations. The step size defines how many integration steps need to be implemented before the final solution is reached. The smaller the step size, the more accurate the solution will be, noting that this will also increase the computing time. The step size corresponds to the characteristic Δh that according to Bellot Rubio et al. [32] can be chosen around a few hundred of meters, 100–300 m.Once the step size is defined, we define the model coefficients as follows:(16)l1=ftn,yn,l2=ftn+12p,yn+p12l1,l3=ftn+12p,yn+p12l2,l4=ftn+p,yn+pl3.For our case, the function to be studied is(17)fhn,vn=K′ρhe−1/2σvo2−v2−1/3⋅v⋅1cosz,and the coefficients are calculated as(18)l1=fhn,vn,l2=fhn+12p,vn+p12l1,l3=fhn+12p,vn+p12l2,l4=fhn+p,vn+pl3.Once the coefficients (equation (16)) are calculated, we compute the solution for the point yn+1 using the following formula:(19)yn+1=yn+16pl1+2l2+2l3+l4.For our case (equation (18)), this is(20)vn+1=vn+16pl1+2l2+2l3+l4.The result obtained is the solution for the point (hn+1, vn+1), which becomes the initial condition to find the numerical approximation for the next point. The procedure is repeated until the desired value is reached.


## 3.1.2. Validating the Code Using the Jacchia Catalog

Once the procedure is defined, we need to validate the code by comparing the results with the previously published data. We use a catalog of very precise photographic trajectories of meteors [58], hereafter referred to as the JVB catalog, obtained using high spatial resolution Super-Schmidt cameras. Jacchia’ [58] study synthesizes the physical inferred parameters for 413 meteors between −5 and +2, 5 stellar magnitude ranges obtained during a multistation meteor network operated in the fifties and sixties in New Mexico, USA. The data provide the meteor velocity and magnitude as a function of the altitude, the derived preatmospheric velocity, the deceleration, and some additional information for the observed meteors. All JVB catalog events were named using a number. In this project, we have used the same Jacchia and Whipple [58] numbering, but included a J at the beginning for the computational reasons. For example, we later discuss meteor J8945 (listed in JVB as 8945).Equation (13) also requires the knowledge of air density. We adopted a general model widely used in meteor studies [45], the United States (US) standard atmosphere [62]. The US standard atmosphere was originally developed in 1958 by the US Committee on Extension to the Standard Atmosphere and improved in 1976. It is a series of tables which approximate the values for atmospheric temperature, density, pressure, and other properties over a wide range of altitudes.


## 3.1.3. Autofitting Procedure

We have defined a way to transform the differential equation into an expression that can be iteratively computed. The aim is to find a result for whichK′ and σ produce the closest curve fitting to the observational data points. In order to find the best matching values, we introduce the following autofitting procedure. We start by choosing two random values of K′ and σ set as initial approximation. The meteor data to be investigated include the velocity of a meteor at specified points of altitude; we compare this velocity with the velocity simulated by our code at the same altitude points. The error factor is then calculated as follows:(21)ξ=vd−vfK,σ2,where vd is the velocity inferred from the measured data, and vf is the computed velocity.Furthermore, we introduce the increment factors forK′ and σ, which are defined as ΔK and Δσ. Following the same procedure as earlier, we compute the error factor for(22)ξ=vd−vfK+ΔK,σ2ξ=vd−vfK−ΔK,σ2,ξ=vd−vfK,σ+Δσ2ξ=vd−vfK,σ−Δσ2,ξ=vd−vfK+ΔK,σ+Δσ2ξ=vd−vfK−ΔK,σ−Δσ2,ξ=vd−vfK+ΔK,σ−Δσ2ξ=vd−vfK−ΔK,σ+Δσ2.In principle, we created a 2D matrix of errors. The error parameter can be shown in a table for better visualization of the algorithm (Figure1). Once all the error values are computed, the search for the minimum value is performed, and the result is considered the new centered value for the next computation iteration.Figure 1
Basic scheme of the error matrix.We repeat the same procedure for the next centered value, until we reach the point where the minimum will be centered in the middle of the matrix. Consequently, the minimum error value will correspond to the sought values ofK′ and σ. If very small increments of K′ and σ are used, the solution will be more accurate (at the cost of increased computing time). If we set large increments of K′ and σ, the solution will be reached faster, but at the expense of the resolution of the solution. In order to deal with this conundrum, the code is optimized such that it is set to work with large increments of K′ and σ at the beginning. Once the ‘first approximation’ solution is found, the code switches to smaller increments until the optimal resolution is reached.Figure2 shows an example of the autofitting procedure. We have chosen the J8945 meteor for the comparison with published data [32]. Other cases were studied and compared in Dergham [29]. By comparing the velocity vs. altitude graphs, it is evident that the fitting curves are very similar to the observed data.Figure 2
Normalized velocity vs. altitude fit for the meteor J8945. The result is a good fit to the previously published Figure1 of Bellot Rubio et al. [32].


## 3.1.4. Abrupt Disruption Cases

We have presented a model capable of obtaining some parameters for meteoroids. However, as previously mentioned, not all meteoroids can be studied using this particular model because if they undergo fragmentation, the results might be skewed. Bellot Rubio et al. [32] also mentioned that there are a significant number of cases in the JVB catalog that cannot be fitted, likely as a consequence of undergoing abrupt disruptions. Figure 3 shows the results of the velocity curve using our model for the case J4141. The results are also compared with that obtained by Bellot Rubio et al. [32].Figure 3
Comparison of results for the meteoroid J4141. The model cannot produce a well-fitting solution since this meteoroid likely suffered an abrupt disruption. This figure can be compared with that shown in Figure1 of Bellot Rubio et al. [32].Despite the difficulties to obtain well-fitting solutions for some events, it is remarkable that our model is able to identify and produce solutions for the events undergoing quasicontinuous fragmentation. A possible way to study meteoroids with abrupt fragmentation is to focus on different segments of the trajectory that do not include a disruption episode. We will apply that approach to the lower part of the trajectory of the Chelyabinsk bolide to demonstrate that the obtained parameters are consistent.


## 3.2. The Ablation Coefficient

Equation (13) has several unknowns. The ability of a meteoroid to produce light can be linked to its mass loss or the ablation coefficient σ [34, 42]. In principle, the ablation coefficient reflects how rapidly the meteoroid loses its mass as it interacts with the atmosphere. Low values of the ablation coefficient indicate that the object is losing a lesser amount of mass compared to that having a higher value of the ablation coefficient. The ablation coefficient is normally represented in units of s2 km−2 and can also be expressed through the dimensionless mass loss parameter [44, 63] used through the range of meteor physics studies. The value of the ablation coefficient depends of many factors, such as chemical composition, grain size, density, porosity, and body shape, among others. In general, the values of the ablation coefficient range between 0.01 and 0.3 s2 km−2 [25]. In order to exemplify this, we applied different ablation coefficients to a 1 g meteoroid with the preatmospheric velocity of 25 km/s and which starts to decelerate at the altitude of 100 km. The results are shown in Figure 4.Figure 4
Comparison of meteoroids having the same initial conditions (m0 = 1 g, v0 = 25 km/s) but different ablation coefficients.In general, the larger the ablation coefficient is, the faster the body decelerates due to a more rapid mass loss. Consequently, the mass of the body decreases due to ablation; this is described by using the ablation coefficient, and the drag force imposed by the atmosphere has a larger effect. Table1 shows the comparison of our results to several events described in the JVB catalog.Table 1
Comparison of the ablation coefficients for various meteor cases from the JVB catalog.
Meteor IDJacchia (s2·km−2)This study (s2·km−2)J68820.08120.075J69590.03310.0382J72160.05010.079J89450.03540.036J90300.05490.0542J71610.03540.0381


## 3.3. Deceleration, Normalized Instantaneous Mass, and Mass Loss Rate

In this section, the effect of deceleration is examined in more detail. Given that we have the velocity as a function of altitude along the trajectory, we can study deceleration from the behavior of the velocity curve. This is particularly useful as many meteor processing algorithms and detection methods provide velocity values sequentially [15, 64]. At a certain point hi, the code computes an increment of velocity over the increment of distance at the points immediately before and after. This can be expressed as(23)dvdhhi≅ΔvΔh=vh2−vh1h2−h1.Considering all the points with known velocity and altitude along the trajectory as input data, equation (23) can be applied directly. Figure 5(a) shows the velocity curve as a function of altitude for the J8945 meteoroid.Figure 5
Curves of the acceleration (a), mass evolution (b), and relative mass loss rate (c) as a function of altitude for the meteoroid J8945.
(a)(b)(c)The normalized instantaneous mass (m/m0) is the next quantity to be studied. The expression for the normalized instantaneous mass can be derived by rearranging equation (9):(24)mmo=e−1/2σvo2−v2.Equation (24) expresses the normalized instantaneous mass as a function of velocity, whereas the values of velocity are the functions of altitude. Figure 5(b) shows the behavior of normalized mass expressed as a function of altitude for the J8945 event.We define the normalized mass loss rate as the derivative of the relative mass over the altitude. This value is computed as follows:(25)dm/moh=Δm/moΔh=m/moh2−m/moh1h2−h1.Figure5(c) shows the direct application of this approach for the J8945 event.


## 3.4. Application to the Chelyabinsk Superbolide

We apply the Runge–Kutta code developed in this work to the famous Chelyabinsk superbolide. On February 15, 2013, it was predicted that the NEA-2012 DA14, discovered a year prior by the Observatory Astronòmic de Mallorca, would approach the Earth within a close minimum distance of only 27700 km. However, while all the attention was focused on anticipating that encounter, another NEA unexpectedly entered the atmosphere over central Asia on February 15, 2013 at 03:20:33 UTC. The bolide disintegrated in the proximity of the city of Chelyabinsk [1]. The Chelyabinsk bolide reached the brightness magnitude of −28, which is brighter than the moon (Figure 6).Figure 6
A casual photograph of the Chelyabinsk bolide taken by Marat Ahmetvaleev. (a) The abrupt fragmentation at a height of 23 km is clearly distinguishable by the main flare. (b) Pictures of the dust trail by the same photographer reveal multiple paths produced by the fragments. Both pictures are courtesy of the author.
(a)(b)As the days passed and the orbits were calculated, scientists discarded a possible association between the two NEAs, as they presented very different heliocentric orbits. Thanks to video cameras (dash cams) found in the majority of Russian motor vehicles and surveillance cameras placed on buildings, the initial trajectory of the bolide was reconstructed and the orbit was determined [2].After the superbolide sightings, many people uploaded various videos to Internet. Since the geographical location of the recorded videos was known, we reconstructed the bolide trajectory, obtaining the values of velocity as a function of altitude. As shown in Table2, the data were obtained from the analysis of our video compilation, the bolide velocity along the terminal part of its trajectory just after the massive fragmentation event taking place at the altitude of 26 km. The step size was determined by the video frame rate, corresponding to differences in height of 200–150 m. Table 2 shows these data.Table 2
Dynamic data of the Chelyabinsk superbolide.
Height (km)V (km/s)Height (km)v (km/s)18.9814.0415.669.7318.7813.8615.539.4618.5813.6815.399.2018.3813.4915.268.9418.1813.2915.138.6817.9913.0915.018.4217.8012.8814.898.1717.6212.6614.777.9217.4412.4414.667.6717.2612.2214.557.4317.0811.9914.457.1916.9111.7514.346.9616.7411.5114.246.7416.5811.2714.156.5216.4211.0214.066.3116.2610.7613.976.1116.1010.5113.885.9215.9510.2513.805.73Figure7(a) depicts the dynamical data after the disruption had occurred at an altitude of 26 km, and the fit is obtained by our model. The graph shows a quite uniform behavior of the Chelyabinsk superbolide after the main fragmentation event.Figure 7
Velocity evolution (a), mass evolution (b), and mass loss rate (c) as a function of altitude for the Chelyabinsk superbolide.
(a)(b)(c)By studying the dynamic curve, the ablation coefficient can be obtained. The derived value is(26)σ=0.034s2⋅km−2.It is quite remarkable that this value, derived for the lower trajectory, provides a similar ablation coefficient as that for the fireballs at much higher altitudes, even though Chelyabinsk was the deepest penetrating bolide ever documented, still emitting light even as it reached the troposphere. Notably, the atmospheric density in these lower regions is about four orders of magnitude higher.The normalized mass evolution for that lower part of the trajectory is shown in Figure7(b), and the mass loss rate evolution is shown in Figure 7(c). In particular, it is certainly encouraging that the model predicts quite well the ablation behavior of the Chelyabinsk bolide in the lower part of its atmospheric trajectory. The light curve of the Chelyabinsk event was normalized using the US government sensor data at a peak of brightness value of 2.7 we1013 W s·r−1, corresponding to an absolute astronomical magnitude of −28 [1]. According to the NASA JPL Chelyabinsk final report [65], the maximum brightness was achieved at an altitude of 23.3 km. This is consistent with our results as the maximum inferred value of mass loss rate from our dynamic data that occur at an altitude of ∼23.5 km (Figure 7(c)).It is well known that the maximum brightness is achieved shortly after the meteoroid catastrophically breaks apart due to the fragmented/pulverized material being exposed to the heat generated by the resulting shock wave.An interesting conclusion that can be obtained directly from these results is the importance of the atmosphere. As mentioned before, the faster the meteoroid, the more rapid the ablation process. Thus, the atmosphere could effectively shield the Earth from very fast impacts since such objects are preferentially and more quickly ablated. However, the less favorable cases are very large objects (especially if their preatmospheric velocity is low), such as the Tunguska impactor that produced an airburst over Siberia in 1908 when it entered the atmosphere at a velocity of ∼30 km/s [66]. The Chelyabinsk superbolide had the preatmospheric velocity of 19 km/s, and its trajectory was consistent with a grazing geometry (high zenith angle).To exemplify this, we have plotted the entry of the Chelyabinsk superbolide for different initial velocities (Figure8). Of course, disrupting NEOs can produce fragments tens of meters in size. If these fragments encounter our planet under certain geometric conditions, they might become a significant source of damage on the ground and casualties. Thus, identifying the existence of asteroidal complexes in the near-Earth environment is crucial for better assessment of impact hazard.Figure 8
Simulation of the Chelyabinsk superbolide with different initial velocities (v = 17.1, 23, and 27 km/s).If the Chelyabinsk superbolide entered at a higher velocity, it would have slowed down faster due to a more rapid mass loss. According to our model, the maximum brightness of a meteoroid occurs when the mass loss reaches the peak. Consequently, the model predicts that a Chelyabinsk-like asteroid moving at a lower velocity could have more destructive potential on the Earth’s surface (Figure8). On the other hand, a projectile of the same rocky composition and moving at higher velocity could have undergone a much faster ablation process, culminating in an explosive airburst, similar to that occurring over Tunguska. Such conclusions imply that the efficiency of the Earth’s atmosphere to shield us from dangerous tens of meters in size asteroids that strongly depend on the relative velocity of the encounter with our planet.It is important to remark that new improvements in detection of fireballs from space could provide an additional progress in studying the luminous efficiency of bolides [34]. In fact, several detections of the Chelyabinsk bolide from space have been reported [67, 68]. Future studies of common events detected from both the ground and space could constrain the role of the observing geometry in signal loss and possible biases in the subsequent determination of velocity, radiated energy, and determination of the orbital elements.


## 3.5. Implications for Impact Hazard

Considering the destructive nature of large extraterrestrial objects (e.g., [24]), it is crucial to identify the existence of asteroidal complexes in the near-Earth environment. Indeed, the Chelyabinsk scale or larger extraterrestrial bodies (10s of meters in size) are of special interest to the planetary defense community because depending on their orbital geometry and impact angle, these objects could pose significant hazard for humans and infrastructure on the ground. Thus, studying the near-Earth environment along with the well-documented events such as the Chelyabinsk superbolide can shed more light on the dynamic processes, as well as fundamental properties of these objects.From the study of meteorite falls and the relative absence of small impact craters [19], we know that meter-sized asteroids are efficiently fragmented when they penetrate into the stratosphere at hypervelocity. The loading pressure in front of the body produces the rock fracture when it overcomes its tensile strength. As a natural consequence, meteorite falls often deliver tens to hundreds of stones, just after this type of break-ups [69, 70]. We have previously described such behavior when discussing the evolution of Chelyabinsk. Additionally, there is a consensual view that the terrestrial atmosphere behaves as an efficient shield for meter to 10s of meters in diameter projectiles. Despite of this, we still need to increase our statistics as more rarely, and probably for low-velocity projectiles under favorable geometric circumstances, crater excavation is still possible. A good example of this was the so-called Carancas event: an impact crater excavated in the Perú’s Altiplane by a one meter-sized chondritic meteorite [71]. Being a quite unusual impact event, it is also likely that a significant number of these events are rarely studied because they happen in remote locations and remain unnoticed. It is obvious that Chelyabinsk and other well-recorded events (e.g., by fireball networks) provide an opportunity to understand in the required detail the behavior of meter-sized meteoroids and their capacity to directly cause human injuries and even casualties.


## 4. Conclusions

We have developed a numerical model employing the Runge–Kutta method to predict the dynamic behavior of meteoroids penetrating the Earth’s atmosphere based on the meteor physics equations [39]. To test the numerical model, we successfully applied it to several meteor events described in the scientific literature. Having validated the numerical model, we studied the deceleration behavior of the Chelyabinsk superbolide in the lower part of its atmospheric trajectory, just at the region following the main fragmentation event where such an approach is applicable. This scheme represents a novel way to study complicated meteoric events by examining only the portion of the trajectory during which the object does not undergo abrupt fragmentation.Our study of the deceleration profile of the Chelyabinsk superbolide has allowed us to reach the following conclusions:(a)
Our numerical model, successfully applied to the lower part of the fireball trajectory, predicts well the main observed characteristics of the Chelyabinsk superbolide. This is quite remarkable as the lower trajectory studied here has a similar ablation behaviour compared to fireballs at higher altitudes. It should be noted that the Chelyabinsk event is the deepest penetrating bolide ever documented, borderline emitting light as it reached the troposphere. Thus, our approach offers a promising venue in studying complex meteoric events in a streamlined and simplified manner.(b)
The best fit to the deceleration pattern provides an average ablation coefficientσ = 0.034 s2·km−2, which is in the range of those derived in the scientific literature(c)
The ablation coefficient is considered constant within each studied trajectory interval. This simplistic approach is probably one of the reasons why this model is not applicable to the entire trajectory of the meteoroids suffering significant fragmentation and catastrophic disruptions. In any case, for the cases studied here, the ablation coefficients obtained in our work are consistent with those reported in scientific literature.(d)
A comparison of the fireball main parameters inferred from both the ground and space could constrain the role of the observing geometry in signal loss and biases in the determination of the velocity, radiated energy, and the orbital elements(e)
NEOs disrupting in the near-Earth region can produce fragments tens of meters in size, which, if encountering our planet under the right geometric circumstances, might be a significant source of hazard for humans and infrastructure on the ground. Consequently, we suggest that identifying the existence of asteroidal complexes in the near-Earth environment is crucial for a better assessment of impact hazard.(f)
Finally, we should not underestimate the hazardous potential of small asteroids as our model indicates that the ability of the Earth’s atmosphere to shield us from such objects strongly depends on the relative velocity of the encounter with our planet



---

*Source: 8852772-2021-03-27.xml*