© ESO, 2017

1. Introduction

The dust content of air is able to strongly influence the state of the high atmosphere and the processes which occur there. A number of phenomena that take place in the high atmosphere are associated with cosmic dust influx, such as an increase in intensity of the scattered solar radiation field, natural radiation and resonance scattering in the lines of different atoms and ions, and the specular reflection of radio waves from ionized meteor trails (Lebedinets 1981; Bronshten 1983). Cosmic dust particles also serve as the condensation nuclei in the noctilucent cloud formation (Bronshten 1984). It is known that heating even a small number of dust particles (up to 5 × 10⁴) can generate an upward air flow with an ascending velocity of 1 cm s^-1; the dust cloud at mid day can create a stratospheric cyclone, as well (Astapovich 1958; Hsu et al. 1997). The influence of aerosol on climate changing predictions is reflected in the review of Kahn (2012).

Fig. 1 Scheme of the flow in case of outstripping ablation (a) and in case of outstripping melting (b). Scheme of flow in disturbed boundary layer (c): solid thin lines indicate melt particles’ streamlines; thick dashed lines indicate undisturbed boundary layer limits. AB is the curvature radius of disturbed streamlines; Fcf is the centrifugal force acting at moving melt particle.

The kinetics of these processes depend strongly on the dust volume concentration, as well as the number distribution of particles by sizes. These quantities are determined by the underlying physics of the particle generation process. Ablation of meteoroids penetrating the Earth’s atmosphere at hypersonic velocities is one of the origins of the dust particles (Bronshten 1983; Ceplecha et al. 1998; Williams 2002). Meteoroids are responsible for the influx of thousands of tons of tiny particles into the Earth’s upper atmosphere annually (Mathews et al. 2001; Janches & Chau 2005). The passage of high-velocity meteoroids through the atmosphere is believed to affect the overall aerothermochemistry.

Ablation is the most complicated process in the meteor phenomenon, in which the mass loss kinetics depends on all the parameters used in calculations (Stober & Jacobi 2008). The simple theory of meteors can evaluate only the total ablating mass influx (Bronshten 1983). However, the most important characteristics of ablation are the particle sizes, the frequency at which the particles are released, number density and size distribution, which influence processes in the meteor wake and the whole meteor phenomenon. These parameters are determined by the mechanics of the particle formation and release from the meteoroid body, so it is of prime importance to investigate the release mechanism and to quantify the particle generation. We investigate meteoroid ablation through modelling constrained by observations, as direct study of the meteoroid fragmentation and the resultant particle size distributions are not practical.

There are several mechanisms that contribute to meteoroid ablation, including evaporation, sputtering, fragmentation, disruption, shelling, spallation, husking (Borovička et al. 2007; Bronshten 1983), and dispersion of the unstable fusion film (Girin 1992). Moreover, the meteoroid fragmentation can proceed in several hypothetical ways (Borovička et al. 2007; Bronshten 1983): progressive fragmentation into parts, which continue to crumble, quasi-continuous detachment of small particles from the main body, and sudden disruption into a number of particles.

It is generally recognized that quasi-continuous fragmentation is one of the main modes of ablation (Lebedinets & Shushkova 1970). Mechanisms of quasi-continuous fragmentation include hydrodynamic instability of the liquid fusion film formed on a solid meteor body (Bronshten 1983), which has not been investigated in detail. Careful consideration of the aerodynamic interaction between an air stream and molten layer of meteoroid has led us to the conclusion that the high-frequency droplet spraying from the molten layer due to the hydrodynamic instability of the melt surface is actual. This hypothesis was originally declared in Girin (1992) and it was based upon the similarity in hydrodynamics of the meteoroid fusion film and a liquid drop in a speedy gas flow. Indeed, values of the hydrodynamic forces, as well as the main criterion for a liquid dispersion (i.e. Weber number), are close for meteoroids and for liquid drops. In its turn, drop shattering is well-studied theoretically and experimentally as a dispersion of the smallest daughter droplets from the parent drop surface (see, for instance, Reinecke & Waldman 1975).

The instability of a thin molten film of a meteoroid was examined in (Girin 1992; Girin & Kopyt 1994) as a mechanism of melt spraying, so that the estimations for the droplet sizes and frequency of their release were obtained. Though, such important characteristics of the meteor phenomenon as transient distributions of the released droplets by sizes, the meteoroid ablation and deceleration laws were not found. Similarly, the ablation of small planetesimals due to the unstable dispersion of fusion film was used as the basic mechanism of spherical chondrules formation by Genge (2000).

The mechanism of the meteoroid melt spraying due to the hydrodynamic gradient instability is considered in the present paper for the particular class of molten meteoroids, leading to the closed mathematical model of the sprayed mass losses. Peculiarities of the gradient instability are described in Sects. 2 and 3; determination of the sprayed particles’ number and sizes are given in Sect. 4, while the model results for three specific meteoroids are considered in Sect. 5. Supplementary material, including the nomenclature list, is in the Appendix.

Generally, the melt instability performance is determined by the balance between meteoroid heating and aerodynamics. Depending on meteoroid velocity, size and substance, two regimes of meteoroid ablation due to instability are possible. When the aerodynamic interaction dominates over the rate of heating, the mass loss overtakes fusion. Conversely, melting outstrips mass losses when the heating dominates over aerodynamics. In the former case the rate of mass loss exceeds the rate of melting, so that the molten layer thickness b is less than the thickness δ_m of the velocity boundary layer in the melt should be. The unstable disturbances are then affected by the stabilizing influence of the meteoroid rigid core (Fig. 1a). In this case, droplets that result from disruption of the molten layer have a size corresponding to a thickness of the layer itself. The results obtained in Girin & Kopyt (1994) are thus valid for the case of dominant ablation.

In the opposite case of the dominant fusion, the molten layer is essentially thicker than the velocity boundary layer induced by the airstream viscous stresses in the melt: b>δ_m. Then, the velocity profile in the melt acquires the form of the completed boundary layer obeying the boundary layer relationships (Schlichting 1960). Thus, the flow near air–melt interface becomes different in this case, consisting of free conjugated boundary layers, in air, δ_a, and in the melt, δ_m, independent from the solid core Fig. 1b). This fact changes principally the mechanism of the instability, compelling the fragmented particles have a size distribution that must be modelled.

2. The instability of a gradient flow near liquid-gas interface

The flow near the liquid − gas interface depends on the media properties. The effect of the velocity profile peculiarities on the liquid − gas interface instability was studied in (Girin 1985). It was established that variations of the media inertial and viscous properties modify the profile, transforming the flow within the conjugated liquid – gas boundary layers. For a linear approximation of the profile (Fig. 2), it was shown that the instability mechanism can differ from the classical Kelvin − Helmholtz type. Namely, the latter type is in action only when the liquid kinematic viscosity is greater than that one of the gas: ν_m>ν_a. In this case the liquid boundary layer is much thicker than the gas one, while the surface velocity, V_s, is much less than the airstream velocity: V_s ≪ V_a. Thus, the velocity profile is close to the discontinuous one (tangential discontinuity) which is the base for the Kelvin − Helmholtz mechanism (Fig. 2a, see Eqs. (2)).

Fig. 2 Linear fit of velocity profile in conjugated boundary layers; case νm>νa(a) and νa>νm(b); δa, δm – air and melt boundary layer thicknesses.

But when the inverse inequality takes place, ν_a>ν_m, the boundary layer thicknesses become comparable, δ_m = O(δ_a), as well as the velocities: V_s = O(V_a) (see Eqs. (2)). Thus, the velocity profile becomes inflated (Fig. 2b) changing the dispersion relation of the boundary-value problem for disturbances. In this case the mechanism of instability becomes completely different from the Kelvin − Helmholtz one, because it is determined now by the unstable root of another dispersion relation, given below as Eq. (1). The instability investigation (Girin 1985) mentioned above resulted in the algebraic dispersion relation of the fourth power with respect to the frequency ω of the disturbance of the form E = E₀exp(ihx − iωt). In the case of liquid − gas systems (i.e. α = ρ_a/ρ_m ≪ 1) it can be simplified to the third-power equation: $\mathrm{(}\mathrm{\Omega }\mathrm{-}\mathrm{\Delta }\mathrm{)}\left[\mathrm{(}\mathrm{\Omega }\mathrm{-}\mathrm{\Delta }\mathrm{)}\mathrm{(}\mathrm{\Omega }\mathrm{+}\mathit{K}\mathrm{\Delta }\mathrm{)}\mathrm{+}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{K}\mathrm{)}\mathrm{\Delta }\right]\mathrm{=}{\mathrm{\Delta }}^{\mathrm{3}}{\mathrm{We}}_{\mathrm{s}}^{-1}\mathrm{(}\mathrm{\Omega }\mathrm{-}\mathit{K}\mathrm{\Delta }\mathrm{)}\mathit{,}$(1)here Ω = ωδ_m/V_s, V_s – velocity at the liquid-gas interface (Fig. 2), ${\mathrm{We}}_{\mathrm{s}}\mathrm{=}{\mathit{\rho }}_{\mathrm{m}}{\mathit{V}}_{\mathrm{s}}^{\mathrm{2}}{\mathit{\delta }}_{\mathrm{m}}\mathit{/}\mathrm{\Sigma }$ is “surface” Weber number, ρ_{a , m} are the air and melt densities, Σ is surface tension of the melt, Δ = 2πδ_m/λ is dimensionless wavenumber, λ is disturbance wavelength, K = [1 − exp( − 2Δ)] / (2Δ), α is the air-to-melt density ratio.

The unstable root of Eq. (1) defines a different type of periodical instability, called the gradient instability (Girin 1985; 2012), the mechanism of which is generated within the liquid boundary layer by the huge velocity gradient, ∂V_{m x}/∂y ≈ V_s/δ_m, which ranges up to (10⁴ − 10⁵)s^-1 for liquid drops as well as for meteoroids. This will be demonstrated later with Eqs. (2). The cause of the gradient instability is the inertial (centrifugal) forces of disturbed curvilinear motion (see Fig. 1c), as explained in Appendix A. As the dispersion relation Eq. (1) contains none of the air parameters ρ_a,V_a,δ_a, the air flow has a weak influence on the disturbances, just transferring momentum to the liquid by means of viscous shear stresses, as it was noted in (Gordillo et al. 2001; Girin 2011a).

The gradient instability mechanism explains the “stripping” mode of a liquid atomization as a quasi-continuous, high-frequency periodic spraying from the unstable part of a liquid surface (Girin 2011a,b). It is known for this mechanism that ablation largely increases the evaporated mass due to the rapid increase of the total liquid surface of the sprayed particles (up to 100 times, Girin 2011b), which would be a significant effect for meteors.

A certain class of middle-sized meteoroids exist, which melt completely because meteoroid heating dominates over ablation. The initial radii of these meteoroids range from 3 × 10^-3 cm to 0.7 cm for iron and from 6 × 10^-3 cm to 0.3 cm for stony meteoroids (Bronshten 1983). Meteoroids larger than these are melted to some depth less than the radius; smaller sized meteoroids decelerate faster than they are totally melted and evaporated (Bronshten 1983). On the other hand, the conclusions of Rogers et al. (2005) show that the smallest of these meteoroids can lose a significant part of their mass due to sputtering. This turns out to have a negligible effect on our meteoroid size considerations, as Fig. 6 of Rogers et al. (2005) shows that sputtering can be neglected for meteoroids with initial radii greater than R ₀> 8 × 10^-3 cm for stony meteoroids, and R ₀> 5 × 10^-3 cm for iron meteoroids. At the same time, the upper size limit for the meteoroids modelled here can also be raised, because the condition b>δ_m does not require that the meteoroid be totally molten. Consequently, in this paper we will consider initial radii ranging from 5 × 10^-3 cm to more than 1 cm for iron and from 8 × 10^-3 cm to 0.5 cm for stony meteoroids.

Thus, the theory of the molten meteoroid spraying can be elaborated in the same way as the theory of drop atomization in a high-speed airflow was constructed (Girin 2011a,b, 2014b). The gradient instability theory is applied below to the class of molten meteoroids in the form of main relationships, which are listed in the next section. A direct extension of this theory, which was validated for a drop in a steady uniform airstream, to meteoroid flight is impossible, as the air density variation with altitude must be taken into account. Therefore, we will consider here the case of a shallow entry angle trajectory (i.e. a small angle γ between the trajectory and the horizon), when the ambient air density is changing gradually with time. In this case, γ ≪ π, we are able to calculate the main kinetic parameters of the sprayed melt droplets: radii r and period t_per of their release from the meteoroid surface, which govern the meteor phenomenon as a whole and the meteor aerosol wake formation in particular.

The meteoroid trajectory in the Earth’s atmosphere can be divided onto two parts: a heating stage, when evaporation and sputtering can be neglected, and an ablation stage, when a meteor itself already exists (Krinov 1960); only the latter stage is considered here. Meteoroid flight proceeds with hypersonic velocities, when the Mach number is large (Ma_∞ ≫ 1), but the air flow behind the detached shock wave is subsonic due to the air compression (Courant & Friedriechs 1976), so that Ma_a ≈ 0.4 (Loitsyanskiy 1995) (Fig. 3). Moreover, the vapour influx increases density in the vicinity of a meteoroid by more than 100 times (Boyd 2000; Popova et al. 2000), further decreasing the Mach number. The role of compressibility is increased only at the lateral sides of a meteoroid, where the gas-stream velocity has some growth. This means that the incompressible medium model is applicable here as a first-order approximation since the compressibility exhibits itself only when the Mach number becomes comparable with unity and greater.

More important is the fact that the gas phase density growth due to the evaporation decreases Knudsen number drastically. Figure 9 of Boyd’s (2000) paper shows that vapour density exceeds the air density near the meteoroid surface by ≈500 times for H-chondrite meteoroid ablating below ≈100 km. That text also states: The mean free path of the flow immediately adjacent to the meteoroid is about 10^-5 m giving a local Knudsen number based on the meteoroid diameter of 0.001 which is in the near-continuum regime. So, following Boyd, the vapours together with the air compression are solely able to diminish the local Knudsen number for meteoroid by 4000 times, which raises the height of the continuum regime by ≈30 km (see Fig. 1 of Campbell-Brown & Koschny 2004 paper). Our usage of the continuum model is thus in agreement with Boyd’s findings. Some other studies, for example, Zinn et al. (2004) use the fluid hydrodynamics equations up to 90 km.

Fig. 3 Flow scheme around meteoroid body. Stable and unstable parts of molten meteoroid surface are shown.

3. The mechanics of particle release

Despite large values of the melt dynamic viscosity, μ_m, the inequality ν_a>ν_m takes place for meteoroids due to the low air density at meteor altitudes. Therefore, the real velocity profile across surface of a completely molten meteoroid is close to that shown in Fig. 2b. For this reason, the gradient instability of the flow in the conjugated air-melt boundary layers is assumed as the mechanism of a particle breakaway from the molten meteoroid surface. The model of a molten meteoroid spraying is built below on the basis of this assumption.

The process of liquid particle generation at some arbitrary elementary area Δl_i is shown in Fig. 4, which is the result of growth of the unstable wave amplitude. It illustrates the natural assumption that the size of a released droplet is determined uniquely by the wavelength λ, while the periodicity of the fragmentation process is given by the characteristic time t_per ≈ Im^-1(ω). When applying theoretical results of an instability treatment to practical problems, the fastest unstable wave is of prime interest, as it is the very wave that will separate a droplet from the meteoroid melt layer. It is sufficient thus for the droplet release description to use the values of the fastest disturbance wavelength, λ_f, and characteristic time t_f (or the growth rate factor, $\mathrm{Im}\mathrm{\left(}{\mathit{\omega }}_{\hspace{0.17em}\hspace{0.17em}\mathrm{f}}\mathrm{\right)}\mathrm{=}{\mathit{t}}_{\mathrm{f}}^{-1}$, Chandrasekhar 1960).

Fig. 4 Particle release process at ith area element of meteoroid surface.

Within the gradient instability theory, these parameters are uniquely determined by the unstable periodic root Ω_f of the dispersion relation Eq. (1). In view of the one-parameter character of the Eq. (1), it appears that the fastest wave parameters only depend on the surface Weber number, We_s. The dimensionless wavenumber Δ_f = 2πδ_m/λ_f and growth rate factor Im(Ω_f) were calculated in (Girin 1985) and are given in Fig. 5 as functions of We_s. The plots show that there exists a critical value We_s.cr. = 3.08 such that the air–melt interface is unstable if We_s> We_s.cr.. At We_s ≥ 20 the parameters become independent from We_s and equal to Δ_f = 1.225, Im(Ω_f) = 0.24. Additionally, the relation between the real and imaginary parts of the unstable root, Re(Ω_f) = 1.5Im(Ω_f), means that periodicity of the particles release is comparable with the characteristic time of e-fold growth of the disturbance amplitude. The behaviour of the wavenumber and growth rate factor depicted in Fig. 5 are applied below to calculate the molten particle release mechanics.

The unstable wave growth depends essentially on the media interaction at the liquid – gas interface. In high-speed flows the interaction occurs within the conjugated boundary layers. On the basis of the general boundary layer theory (Schlichting 1960), the relations between the main boundary-layer flow parameters δ_m, δ_a, V_s, V_a (Fig. 2) can be derived (see Appendix B), which are important in the gradient instability analysis: $\begin{array}{ccc}{\mathit{\delta }}_{\mathrm{m}}\mathrm{=}{\left(\frac{\mathit{\alpha }}{{\mathit{\mu }}^{\mathrm{2}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{3}}{\mathit{\delta }}_{\mathrm{a}}\mathit{,}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathit{V}}_{\mathrm{s}}\mathrm{=}\frac{{\left(\mathit{\alpha \mu }\right)}^{\mathrm{1}\mathit{/}\mathrm{3}}}{\mathit{p}}{\mathit{V}}_{\mathrm{a}}\mathit{,}& & \end{array}$(2)where μ = μ_a/μ_m is air-to-melt viscosity ratio, p = 1 + (αμ)^{1 / 3}.

The flow geometry around the meteoroid is determined by the meteoroid body shape. Assuming the spherical shape of the body and the velocity potential of the air flow around a sphere in the form V_a(ϕ) = 1.5(V_∞ − w)sinϕ (Kochin et al. 1964), we can apply the boundary layer thickness distribution along the meteoroid surface in Ranger’s (1972) form: ${\mathit{\delta }}_{\mathrm{a}}\mathrm{\left(}\mathit{\varphi ,t}\mathrm{\right)}\mathrm{=}\mathrm{2.2}\mathit{R}\mathrm{\left(}\mathit{t}\mathrm{\right)}{\mathrm{Re}}_{\mathrm{a}}^{-0.5}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{\Psi }\mathrm{\left(}\mathit{\varphi }\mathrm{\right)}$. The latter is valid for a gas over a liquid sphere, here Ψ(ϕ) = ((6ϕ − 4sin2ϕ + 0.5sin4ϕ) / sin⁵ϕ)^0.5, ϕ is the polar angle of a surface area element (Figs. 1b and 3). Together with Eq. (B.1) this yields the geometry factor for a liquid sphere: $\mathit{G}\mathrm{\left(}\mathit{\varphi ,t}\mathrm{\right)}\mathrm{=}\mathrm{1.1}\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{R}\end{array}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{\Psi }\mathrm{\left(}\mathit{\varphi }\mathrm{\right)}$, $\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{R}\end{array}\mathrm{=}\mathit{R}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathit{/}{\mathit{R}}_{\mathrm{0}}$. Taking into account Eq. (2), the necessary condition for the gradient instability performance, We_s ≥ 3.08, can be written at an arbitrary area element ϕ of the meteoroid surface in the following form (see Appendix B): ${\mathrm{We}}_{\mathrm{s}}\mathrm{\left(}\mathit{\varphi ,\tau }\mathrm{\right)}\mathrm{=}\frac{\mathrm{2.475}}{{\mathit{p}}^{\mathrm{2}}}\sqrt{\begin{array}{c}\mathrm{˜}\\ \mathit{R}\end{array}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}{\left(\mathrm{1}\mathrm{-}\mathit{W}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\right)}^{\mathrm{3}}}{\sin }^{\mathrm{2}}\mathit{\varphi }\mathrm{\Psi }\mathrm{\left(}\mathit{\varphi }\mathrm{\right)}\hspace{0.17em}\mathrm{GI}\mathrm{\ge }\mathrm{3.08.}$(3)Here $\mathrm{GI}\mathrm{=}{\mathrm{We}}_{\mathrm{\infty }}{\mathrm{Re}}_{\mathrm{\infty }}^{-0.5}$ appears as the gradient instability criterion, ${\mathrm{We}}_{\mathrm{\infty }}\mathrm{=}{\mathit{\rho }}_{\mathrm{\infty }}\mathrm{2}{\mathit{R}}_{\mathrm{0}}{\mathit{V}}_{\mathrm{\infty }}^{\mathrm{2}}\mathit{/}\mathrm{\Sigma }$ and Re_∞ = ρ_∞2R₀V_∞/μ_∞ are the Weber and Reynolds numbers for the meteoroid, w is the meteoroid velocity deficit due to deceleration, so that V_∞ − w is the relative air – meteoroid velocity, W = w/V_∞, τ = t/t_ch; t_ch = 2R₀/α^0.5V_∞ is the characteristic timescale, V_∞ is meteoroid entry velocity.

Applying the function behaviour shown in Fig. 5 to the meteoroid surface, we find out that there exists a critical angle along the surface, ϕ_cr(t), where We_s(ϕ_cr) = We_s.cr. = 3.08. This critical point divides the surface into stable, ϕ<ϕ_cr, and unstable, ϕ>ϕ_cr, parts (Fig. 3). Equality in Eq. (3) defines the value of ϕ_cr(τ). As the normalized meteoroid radius, $\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{R}\end{array}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}$, and the relative meteoroid speed, 1 − W(τ), decrease with time, the initial moment τ = 0 (when $\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{R}\end{array}\mathrm{=}\mathrm{1}\mathit{,W}\mathrm{=}\mathrm{0}$) is the most favourable for the fulfillment of inequality (3). At τ = 0 we will have ϕ_cr<π/ 2 in Eq. (3) when $\begin{array}{ccc}\mathrm{GI}\hspace{0.17em}\hspace{0.17em}\mathit{>}\hspace{0.17em}\hspace{0.17em}\mathrm{0}\mathit{.4}\mathit{.}& & \end{array}$(4)This indicates that the part of the surface adjacent to the meteoroid edge, ϕ_cr<ϕ<π/ 2, is unstable at GI > 0.4, providing a possibility of the high-frequency periodic spraying of the molten substance from the meteoroid surface. This is depicted in Fig. 3. Equation (3) indicates that when GI increases, ϕ_cr decreases, so for large GI ≫ 0.4, the values of ϕ_cr are small, ϕ_cr ≪ π (see the calculated data in Table 1). Thus, the majority of the meteoroid surface generates a mist of the molten particles quasi-continuously, as shown in Fig. 4.

The inequality ${\mathrm{We}}_{\mathrm{\infty }}{\mathrm{Re}}_{\mathrm{\infty }}^{-0.5}⪆\mathrm{0.3}$ was first found empirically by Rabin et al. (1960) as a condition for the stripping mode of a liquid drop breakup. For the standard air conditions, α = 1.3 × 10^-3, it was revealed and grounded theoretically in (Girin 1985) as the criterion for the gradient instability existence at a drop surface. While Eq. (3) expresses the critical conditions locally at every meteoroid surface area ϕ, the criterion GI determines the possibility of instability existence in a whole meteoroid flight.

The inequality Eq. (4) reflects the critical condition for air-melt instability rigorously. We note that the near-critical conditions in the vicinity of the critical point at meteoroid surface yield long wavelengths, which are greater than the meteoroid radius, coupled with large values of the induction period of the instability, which is greater than the period of meteor existence. In practice, this suggests modifying the condition of instability up to GI ≥ 0.6 (or equally We_s ≥ 4.62), see Girin (2011a).

The two Weber numbers, surface We_s and meteoroid We_∞, have the same general sense: the ratio of the hydrodynamic forces to the surface tension forces. But the surface Weber number We_s specifically relates hydrodynamic pressure ${\mathit{\rho }}_{\mathrm{m}}{\mathit{V}}_{\mathrm{s}}^{\mathrm{2}}\mathit{/}\mathrm{2}$ of the melt motion in the liquid boundary layer (which defines the centrifugal force scale) to the surface tension forces δ_m/ Σ of the disturbed interface, whose curvature radius R_c is proportional to δ_m (Fig. 1c). The former force creates the destabilizing effect, defining the centrifugal force scale; the latter force defines a stabilizing effect. Thus, the surface Weber number reflects locally, at a disturbance level, the ratio of destabilizing forces to the stabilizing ones, governing the disturbance development.

The existence of the critical value We_s.cr. has a simple physical meaning, which is explained in Appendix C. This critical value, when applied to a meteoroid, yields the critical hydrodynamic pressure of the airstream running against the meteoroid, $\mathit{P}\mathrm{=}\mathrm{0.5}{\mathit{\rho }}_{\mathrm{a}}{\mathit{V}}_{\mathrm{a}}^{\mathrm{2}}$, which is sufficient to pull out the droplets from the meteoroid melt by means of the unstable disturbances. These values are calculated in Appendix C for the three particular meteoroids considered in the present paper: ${\mathit{P}}_{\mathrm{cr}}^{\mathrm{I}}\mathrm{=}\mathrm{3.0}\mathrm{kPa}$, ${\mathit{P}}_{\mathrm{cr}}^{\mathrm{II}}\mathrm{=}\mathrm{2.1}\mathrm{kPa}$, ${\mathit{P}}_{\mathrm{cr}}^{\mathrm{III}}\mathrm{=}\mathrm{0.9}\mathrm{kPa}$.

Fig. 5 Basic dependencies Δf(Wes) (solid line) and Im [ Ωf(Wes) ] (dashed line) of gradient instability.

Previously, it has been noted (Borovička et al. 2007) that dynamic pressure of the airstream does not exceed the meteoroid strength and therefore the fragmentation cannot be explained by the direct action of the airstream. The dynamic pressure is on the order of 5 − 20 kPa at heights just below 100 km, which is too small to prevail over the meteoroid strength. Conversely, the critical values of $\mathrm{0.5}{\mathit{\rho }}_{\mathrm{a}}{\mathit{V}}_{\mathrm{a}}^{\mathrm{2}}$, derived above, which correspond to the value We_s.cr., show that the dynamic pressure of the airstream at heights up to 100 km, is sufficient to overcome surface tension forces and pull out liquid particles from the meteoroid melt by means of unstable disturbances. Hence, this model is able to explain meteoroid fragmentation at large heights and low dynamic pressures.

On the other hand, the melt spraying mechanism resembles models of quasi-continuous meteoroid fragmentation (such as dustball, thermal erosion, and thermal disruption models, see Campbell-Brown et al. 2013), but is able to provide theoretical predictions of the fragmentation frequency, and fragment size distributions. These are based on the rigorous mathematical treatment of the instability problem, which was carried out in Girin (1985).

4. The quantification of the melt spraying

Because of the small values of ρ_∞ at higher altitudes, Eq. (4) is satisfied only below the altitude z₀, where the air density becomes large enough to produce a mechanical effect. Above that, the meteoroid interaction with the atmosphere has thermal nature and leads to the meteoroid heating and the molten layer formation (Krinov 1960). Values of z₀ are calculated for specific meteoroids and given in Table 1.

The values of z₀ reflect the critical condition in Eq. (4) for gradient instability, but there are some additional necessary conditions. For instance, the meteoroid must first become molten for the gradient instability itself to work. According to Bronshten’s (1983) data for R₀ = 1 cm, the heights z_m, where the meteoroid attains a melt temperature T_m, for the considered in present paper three variants are as follows: ${\mathit{z}}_{\mathrm{m}}^{\mathrm{I}}\mathrm{=}\mathrm{105}\mathrm{km}$, ${\mathit{z}}_{\mathrm{m}}^{\mathrm{II}}\mathrm{=}\mathrm{88}\mathrm{km}$, ${\mathit{z}}_{\mathrm{m}}^{\mathrm{III}}\mathrm{=}\mathrm{112}\mathrm{km}$. In each case we have z_m ≤ z₀, so, the spraying process starts well below z₀, where the melt layer grows sufficiently and the gradient instability mechanism is able to start. The altitude, where the meteoroid is molten through the characteristic depth of its half-radius, can be estimated with the use of the Lebedinets & Portniagin (1966) solution for the temperature distribution inside a rapidly and randomly rotating sphere. According to this solution, the ratio of densities at the heights where the meteoroid surface attains the melt temperature first, ρ_m, and where the meteoroid body is melted by the depth of half radius, ρ_h.r., is: ${\mathit{\rho }}_{\mathrm{m}}\mathrm{=}{\mathit{\rho }}_{\mathrm{h}\mathit{.}\mathrm{r}\mathit{.}}\frac{{\mathit{R}}_{\mathrm{0}}}{\mathrm{0.5}{\mathit{R}}_{\mathrm{0}}}\frac{\sinh \mathrm{(}\mathrm{0.5}{\mathit{R}}_{\mathrm{0}}\mathit{/}{\mathit{x}}_{\mathrm{0}}\mathrm{)}}{\sinh \mathrm{(}{\mathit{R}}_{\mathrm{0}}\mathit{/}{\mathit{x}}_{\mathrm{0}}\mathrm{)}}$, where ${\mathit{x}}_{\mathrm{0}}\mathrm{=}\mathit{b}\mathrm{(}\frac{{\mathit{H}}^{\mathrm{\ast }}}{{\mathit{V}}_{\mathrm{\infty }}\sin \mathit{\gamma }}{\mathrm{)}}^{\mathrm{0.5}}$ is the characteristic depth of heating. So, the ratio ρ_h.r./ρ_m depends on the meteoroid substance, size and velocity. With the use of data for x₀ (Bronshten 1983), we have R₀/x₀ = 3.33;1.76;6.58 consecutively for the meteoroids considered in the present paper. Making calculations with the formula above, we will have: ρ_h.r./ρ_m = 2.75;1.41;6.74 consecutively. Then the difference in the two heights may be determined assuming the uniform atmosphere density distribution with altitude. In this way we obtain: z_h.r. = z_m − 7.1 km;2.4 km;13.4 km consecutively. The height at which the melt depth is sufficient for the spraying process beginning can be estimated as z_h.r. = 97.9 km;85.6 km;98.6 km consecutively for the meteoroids considered in the present paper. Thus, the meteoroid spraying will start at heights smaller than z = z₀ in all of the variants considered here.

Once the critical condition in Eq. (4) is fulfilled for a meteoroid, each elementary area, where the inequality Eq. (3) is satisfied, can be considered as a high-frequency source of the molten meteoroid particles. In our analysis we must consider the air flow variability along the curved meteoroid surface, which causes the dependence of λ _f, t_f and We_s on the polar angle ϕ (Fig. 3). We divide the meteoroid surface into a system of plane elementary areas (more exactly − spherical belts in the axially symmetrical flow) of the width Δl_i = R(t) Δϕ_i (Fig. 4), so that at every small time interval Δt the flow parameters may be treated as constant at every ith area element. The problem of the breakaway droplet parameters determination is then reduced to the application of the above mentioned instability regularities to this local-plane flow. Only the windward part of a meteoroid surface is considered further since the flow past the meteoroid leeward part is weak.

As the droplets are generated by the fastest unstable disturbance, it is natural to assume that the droplet sizes and period of their release from the main body are proportional to the disturbance wavelength and characteristic time of e-fold growth of the disturbance amplitude at a given ith area element, correspondingly: $\begin{array}{ccc}{\mathit{r}}_{\mathit{i}}\mathrm{=}{\mathit{k}}_{\mathrm{r}}{\mathit{\lambda }}_{\mathrm{f}\mathit{i}}\mathrm{=}{\mathit{k}}_{\mathrm{r}}\mathrm{2}\mathit{\pi }{\mathit{\delta }}_{\mathrm{m}}\mathit{/}{\mathrm{\Delta }}_{\mathrm{f}}\mathit{,}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathit{k}}_{\mathrm{r}}\mathrm{\le }\mathrm{0.25}\mathrm{;}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}& & \\ {\mathit{t}}_{\mathrm{per}\mathit{i}}\mathrm{=}{\mathit{k}}_{\mathrm{t}}\hspace{0.17em}{\mathit{t}}_{\mathrm{f}\mathit{i}}\mathrm{=}{\mathit{k}}_{\mathrm{t}}{\mathrm{Im}}^{-1}\mathrm{(}{\mathrm{\Omega }}_{\mathrm{f}}{\mathit{V}}_{\mathrm{s}}\mathit{/}{\mathit{\delta }}_{\mathrm{m}}\mathrm{)}\mathit{,}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathit{k}}_{\mathrm{t}}\mathrm{\approx }\mathrm{1.}& & \end{array}$The model coefficients, k_r = 0.17 and k_t = 1.1, were determined in Girin (2011a) by constraining the model with experimental data.

Having the main parameters of the spraying process, r_i and t_peri, expressed via the main parameters of the air flow, δ_a and V_a in Eqs. (5) and (6) with the use of Eq. (2), we are able to quantify the spraying process. The number of wavelengths which are confined within ith area element, is equal to n_i = Δl_i/λ _{f i}. Due to the axial symmetry of the flow past the meteoroid, this is the number of tori of radius R(t)sinϕ_i, which have been broken away from the spherical belt corresponding to Δl_i. The toroidal breakaway structures were observed in LIF-experiments (Laser Induced Fluorescence) of Theofanous & Li (2008) on a liquid drop breakup.

When the disturbance has grown in amplitude by e, or τ_indi = ^∫dt/t_peri ≥ 1, we suppose that separation of n_tori = n_i tori of cross-section radius r_i occurs at that moment from that area element. Assuming that the liquid tori rapidly break up into droplets of radius r_i in the airstream, and taking the overall volume of the tori, $\mathrm{\Delta }{\mathit{\upsilon }}_{\mathit{i}}\mathrm{=}{\mathit{n}}_{\mathrm{tor}\mathit{i}}\mathrm{·}\mathrm{\Delta }{\mathit{\upsilon }}_{\mathrm{tor}\mathit{i}}\mathrm{=}\mathit{\pi }{\mathit{k}}_{\mathrm{r}}^{\mathrm{2}}{\mathit{\lambda }}_{\mathrm{f}\hspace{0.17em}\mathit{i}}\mathrm{2}\mathit{\pi R}\sin {\mathit{\varphi }}_{\mathit{i}}\mathrm{·}\mathrm{\Delta }{\mathit{l}}_{\mathit{i}}$, divided by the single droplet volume, υ_i = 4 / 3πr_i³, we obtain the differential equation for the number Δn_i of the newborn particles: $\begin{array}{ccc}\mathrm{\Delta }{\mathit{n}}_{\mathit{i}}\mathrm{\left(}{\mathit{\varphi }}_{\hspace{0.17em}\mathit{i}}\mathit{,\tau }\mathrm{\right)}\mathrm{=}{\mathit{B}}_{\mathrm{2}}\sqrt{\begin{array}{c}\mathrm{˜}\\ \mathit{R}\end{array}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}{\left(\mathrm{1}\mathrm{-}\mathit{W}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\right)}^{\mathrm{5}}}\frac{{\sin }^{\mathrm{2}}{\mathit{\varphi }}_{\hspace{0.17em}\mathit{i}}}{{\mathrm{\Psi }}^{\mathrm{3}}\mathrm{\left(}{\mathit{\varphi }}_{\hspace{0.17em}\mathit{i}}\mathrm{\right)}}\mathrm{\Delta }{\mathit{\varphi }}_{\hspace{0.17em}\mathit{i}}\mathrm{\Delta }\mathit{\tau ,}& & \end{array}$(7)which are released during period Δτ = t_indi/t_ch from the area element Δl_i = RΔϕ_i. In accord with Eq. (5) and expression Δ_f = 2πδ_m/λ_f, the newborn particle radius can be written as follows: $\begin{array}{ccc}{\mathit{r}}_{\mathit{i}}\mathrm{\left(}{\mathit{\varphi }}_{\hspace{0.17em}\mathit{i}}\mathit{,\tau }\mathrm{\right)}\mathrm{=}{\mathit{R}}_{\mathrm{0}}{\mathit{B}}_{\mathrm{1}}\hspace{0.17em}\sqrt{\frac{\begin{array}{c}\mathrm{˜}\\ \mathit{R}\end{array}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}}{\mathrm{1}\mathrm{-}\mathit{W}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}}}\mathrm{\Psi }\mathrm{\left(}{\mathit{\varphi }}_{\hspace{0.17em}\mathit{i}}\mathrm{\right)}\mathit{.}& & \end{array}$(8)$\begin{array}{ccc}\mathrm{Here}\hspace{0.17em}{\mathit{B}}_{\mathrm{1}}\mathrm{=}\frac{\mathrm{4.4}\mathit{\pi }\hspace{0.17em}{\mathit{k}}_{\mathrm{r}}}{{\mathrm{\Delta }}_{\mathrm{f}}\mathrm{(}\mathrm{2}{\mathrm{Re}}_{\mathrm{\infty }}{\mathrm{)}}^{\mathrm{0.5}}}\frac{{\mathit{\alpha }}^{\mathrm{1}\mathit{/}\mathrm{3}}}{{\mathit{\mu }}^{\mathrm{2}\mathit{/}\mathrm{3}}}\mathit{,}\hspace{0.17em}{\mathit{B}}_{\mathrm{2}}\mathrm{=}\frac{\mathrm{0.21}{\mathrm{\Delta }}_{\mathrm{f}}^{\mathrm{2}}\mathrm{Im}\mathrm{\left(}{\mathrm{\Omega }}_{\mathrm{f}}\mathrm{\right)}\mathrm{(}\mathrm{2}{\mathrm{Re}}_{\mathrm{\infty }}^{\mathrm{3}}{\mathrm{)}}^{\mathrm{0.5}}}{\mathit{\pi }\hspace{0.17em}{\mathit{k}}_{\mathrm{r}}{\mathit{k}}_{\mathrm{t}}\mathrm{(}\mathrm{1}\mathrm{+}\mathrm{(}\mathit{\alpha \mu }{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{3}}\mathrm{)}}\frac{{\mathit{\mu }}^{\mathrm{7}\mathit{/}\mathrm{3}}}{{\mathit{\alpha }}^{\mathrm{7}\mathit{/}\mathrm{6}}}& & \end{array}$have the sense of scaling parameters for the particle sizes and number, respectively (Girin 2011a). For large values of We_s, We_s> 20, B₁ and B₂ are constants. The derivation of Eqs. (7) and (8) can be found in Appendix D.

Fig. 6 Distributions of breakaway droplets by sizes: Δn(r)(a) and ΔM(r)(b). Variant I: fast iron meteoroid. Blue indicates the final distribution at τ = 2.72; brown, red, light green, and crimson indicate, respectively, the intermediate distributions at τ = 0.44, 0.89, 1.39, and  1.80; tch = 2.16 × 10-3 s. See Movie 1.

Fig. 7 Distributions of breakaway droplets by sizes: Δn(r)(a) and ΔM(r)(b). Variant II: slow iron meteoroid. Blue indicates the final distribution at τ = 3.11; brown, red, light green, and crimson indicate, respectively, the intermediate distributions at τ = 0.30, 0.93, 1.88, and 2.48; tch = 5.15 × 10-3 s. See Movie 2.

In general, the intricacy of Δ_f and Im(Ω_f) dependencies on ϕ and t via surface Weber number We_s(ϕ,t) together with the transiency and complexity of the air flow geometry makes it necessary to carry out numerical study of the ablation process (though in part two of this paper we will present the approximate analytical solution to the problem). As previously done for a drop atomizing in a speedy gas stream (Girin 2014b), the values for V_a,V_s ,δ_a,δ_m, We_s, Δ_f, Im(Ω_f),r,τ_ind and Δn are computed here with Eqs. (2), (3), (5)–(8) and results of Fig. 5. Then, the melt mass $\mathrm{\Delta }\mathit{m}\mathrm{=}\sum _{\mathit{i}}\mathrm{\Delta }{\mathit{n}}_{\mathit{i}}\mathrm{·}\mathrm{4}\mathit{/}\mathrm{3}\mathit{\pi }{\mathit{\rho }}_{\mathrm{m}}{\mathit{r}}_{\mathit{i}}^{\mathrm{3}}$, which is released as breakaway particles, is subtracted from the meteoroid mass m(t) as the outcome of ablation; the current meteoroid radius R(t) is thus found.

Fig. 8 Distributions of breakaway droplets by sizes: Δn(r)(a) and ΔM(r)(b). Variant III: stony meteoroid. Blue indicates the final distribution at τ = 105; purple, orange, light-blue, green, brown, red, light green, and crimson indicate, respectively, the intermediate ones at τ = 2.05, 3.88, 5.94, 9.96, 16.8, 29.9, 48.7, and 73.0;tch = 1.45 × 10-3 s.

5. Results and discussion

The model of melt spraying presented here allows us to determine the most important parameters of the ablation process, including the distribution of all sprayed particles by sizes, as well as evolution of the sprayed mist in space and time. Here we consider three particular cases of spherical meteoroids of initial radius R₀ = 0.3 cm and different velocities and substances, which begin to ablate when moving in the continuum regime. The six-fold air compression and flow deceleration behind the detached shock front ahead of the meteoroid (Fig. 3) are taken into account. The density of air flow and meteoric material around the meteoroid is assumed to be large enough, so that the gas species density value, ρ_a = 10^-7 g cm^-3, is much greater than the ambient atmosphere density ρ_∞ and corresponds to almost 30 km lower height (see the end of Sect. 2).

Results of the calculations in the form of transient number, Δn(r), and mass, ΔM(r), distributions of the sprayed droplets by sizes are given in Figs. 6–8. In addition, two video files, Movie 1 and Movie 2 illustrate the transient dynamics of the particle size distribution for the considered meteoroids. During flight, the meteoroid velocity is assumed to change in accordance with empirical law: W = 1 − exp( − Cτ) (Girin 2011a), where parameter C = 2α^0.5 has the sense of meteoroid deceleration, as the derivative dW/ dτ shows. In this case, the meteoroid velocity decrease is negligible along most of the trajectory till the end of meteoroid ablation, because the characteristic time of meteoroid deceleration is much greater than the one of spraying: C^-1 ≫ τ_ch. The calculated dependence of meteoroid radius on the count-down time, t_s.d. − t, is close to linear, so that the ablation law is close to the form M = (1 − t/t_s.d.)³; here t_s.d. is the overall time of the spraying process duration.

Some general information is given in Table 1, where the initial meteoroid velocity, V_∞, is given in km s^-1 and the ablation coefficient, $\mathit{\sigma }\mathrm{=}{\mathrm{\left(}\mathit{S}{\mathit{\rho }}_{\mathrm{\infty }}{\mathit{V}}_{\mathrm{\infty }}^{\mathrm{3}}\mathrm{\Gamma }\mathrm{\right)}}^{-1}\hspace{0.17em}\mathrm{d}\mathit{m}\mathit{/}\mathrm{d}\mathit{t}$, is in s² cm^-2; here Γ is the drag coefficient. The media properties are: ρ_ir = 7.8 g cm^-3, ρ_st = 3.5 g cm^-3, μ_∞ = 6.8 × 10^-5 kg m^-1 s^-1, μ_ir = 5.8 × 10^-3 kg m^-1 s^-1, μ_st = 0.174 kg m^-1 s^-1, Σ_ir = 1.2 N m^-1, and Σ_st = 0.36 N m^-1. The values of Reynolds criterion Re_∞ indicate that the regime of air streamlining past meteor body provides in each case the development of boundary layer in a liquid melt, while GI-criterion values are large enough to provide a full-strength working of the gradient instability mechanism. The minimum value of the growth time t_f of the disturbance amplitude, taken along the meteoroid surface, is shown as the characteristic time of the breaking-away periodicity of the molten droplets. Compared with the overall spraying time t_s.d., the calculated values of t_f indicate the intensity of the spraying mechanism, since the ratio t_f/t_s.d. is small, ranging from 0.001 to 0.002 in the considered cases of meteoroids.

One can see that the ablation coefficient for the iron meteoroids falls into the range of the values generally accepted in the physical theory of meteors, σ = (10^-12 − 10^-11) s² cm^-2 (Bronshten 1983), while the value for the stone meteoroid is two orders lower than this range.

But this low value correlates with the intrinsic values of ablation coefficient, which were revealed and obtained recently (Ceplecha & ReVelle 2005; Borovička 2006). Indeed, when the fragmentation mode of ablation proceeds, the fragments ablate simultaneously, so that two different values of ablation coefficient, intrinsic and apparent, must be considered. Specifically, the former lies in the range from 10^-13 to 8 × 10^-13 s² cm^-2, while the latter has considerably greater values. The intrinsic ablation coefficient is calculated here.

Calculations reveal a large number, N = 1.5 × 10⁶, of the particles produced by the fast iron meteor, which are two orders smaller by sizes than the parent meteoroid. Their sizes range from 4 μm to 135 μm, with the median radius r_med = 26.9 μm (Fig. 6); however, the mass fraction of those particles which have radii larger than r_med, constitutes 0.76m₀. A few particles have sizes greater than 45 μm. These values can be compared with those obtained in Campbell-Brown et al. (2013) for meteor 20101106_065820 (velocity V_∞ = 69.5 km s^-1, radius R₀ = 0.4 mm) with the use of the thermal disruption model. Specifically, the particle sizes were ranged from 31 μm to 144 μm.

As the expressions for the scaling parameters B₁, B₂ show, the particle sizes vary generally as ${\mathit{V}}_{\mathrm{\infty }}^{-0.5}$, while their total number N varies as ${\mathit{V}}_{\mathrm{\infty }}^{\mathrm{1.5}}$. For the slower iron meteoroid, V_∞ = 25 km s^-1, the sizes are larger and range from 16 μm to 265 μm, with the median radius r_med = 41.6 μm; a few particles have sizes greater than 70 μm. The mass fraction of those particles which have a radius greater than r_med, constitutes 0.78m₀. The dynamics of the transient distributions of the slow and fast iron meteoroids are reflected in the supplementary video files (Movie 1 and Movie 2).

The rate of ablation varies along the meteoroid surface. In all calculations it had its minimum at the angle ϕ ₁ near the critical angle and then increased towards ϕ = π/ 2. The larger droplets are always generated at the beginning of the ablation process, while the smaller droplets are released at the end, when the meteoroid radius is small. This means that, when the molten meteoroid moves along a gently sloped trajectory, the beginning of its aerosol trail consists of more coarse particles than the end.

The melt of iron meteoroid has boundary layer nearly 13 times thinner than that of a stony meteoroid, along with an interface velocity V_s that is 2.5 times larger as formulae in Eq. (2) show. This yields a more that 30-fold increase in velocity gradient in the iron melt (as it is directly proportional to μ, see Eq. (2)) and a consequently stronger destabilization of the iron-melt flow. This consideration leads to the conclusion regarding the significant role of thermal influx to the melt, which must be taken into account in instability applications, as the melt viscosity depends considerably on temperature.

The stony meteoroid disintegrates much more slowly than the iron ones due to the 30-fold larger viscosity, which agrees with conclusions of Bronshten (1983). This leads to large increase in the total ablation time, t_s.d., and to a difference in the breakaway particle sizes of more than an order, in comparison with the iron meteor of the same velocity. The particle sizes range from 100 μm to 570 μm, with the median radius r_med = 336 μm, as seen in Fig. 8. The mass fraction of those particles, which have radii greater than r_med, constitutes 0.74m₀. These values can be compared with those obtained with the use of the thermal disruption model in Campbell-Brown et al. (2013) for meteor 20101103_053624, which had a little greater velocity (V_∞ = 66.4 km s^-1) and less radius (R₀ = 0.7 mm). Namely, the particle sizes ranged from 81 μm to 272 μm.

The fragmentation of the stony meteoroid is not terminated by the gradient instability mechanism. As seen in Eq. (2), this is due to the thick boundary layer, which for a high-viscosity stony melt becomes comparable with the radius of the meteoroid remnant. Thus, the fastest unstable disturbance wavelength, λ_f, becomes comparable with the remnant radius $\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{R}}_{\mathrm{rem}}\end{array}\mathrm{=}\mathrm{0.14}$. We note that this terminal fragmentation has little effect on the overall meteoroid ablation since the remnant meteoroid mass is m_rem = 2.7 × 10^-3m₀. However, the molten remnant can be fragmented afterwards by the mechanism of the Rayleigh − Taylor instability, since the current Weber number value is sufficiently large at that moment, We_∞ ≈ 100, to provide the “claviform” or “multibag” mode of the remnant breakup (Reinecke & Waldman 1975; Zhao et al. 2010; Girin 2012, 2014a). Some ending meteor flashes can be explained thus as the delayed bursting of their remnants (Girin 1992; Girin & Kopyt 1994), since the low-Weber modes of breakup consist in a sudden rupture of a vast thin “bag” film into the tiny droplets.

Two−mode particle distribution may exist at small values of the gradient instability criterion GI due to the influence of the “hump” in Δ_f(We_s) dependence shown in Fig. 5. This feature causes a shift towards smaller particle sizes near the critical point ϕ_cr on the meteoroid surface. At smaller GI values this effect becomes more distinct, as Fig. 7 shows, since ϕ_cr becomes greater at smaller GI (see Table 1), approaching ϕ = π/ 2. As a result, most part of the meteoroid surface produces smaller particles, creating the bi−modal distribution in droplet sizes. The two−mode distribution can be seen in paper Campbell-Brown et al. (2013) in the analysis of low-speed meteor 20101016_070052.

Table 1 shows that the meteor path, l, for the iron meteoroid is shorter than for the stony one. Together with the difference in the total number of generated particles this gives significantly different aerosol density in the meteor wakes.

6. Conclusions

We have presented the theory of a molten meteoroid spraying, which allowed an investigation of the details of ablation kinetics. The theory is based on the mechanism of the gradient instability, which generated high-frequency quasi-continuous spraying of the liquid melt particles constituting a meteor aerosol wake.

The analysis shows that the dynamic pressure of the airstream at heights 80–100 km is sufficient to pull out the liquid melt droplets via the mechanism of gradient instability, which provides a simple and natural explanation of the quasi-continuous fragmentation of small (≈1 cm) meteoroids at similar heights.

The problem of meteoroid ablation is solved theoretically in a closed form for the case of a shallow entry angle, the physical properties, initial radius and velocity of the meteoroid being given. The sizes of molten, sprayed droplets and transient distributions for the iron and stony meteoroids are calculated. Differences in iron and stony meteoroid mechanical properties cause significant changes in the generated particles: meteor wake lengths, particle size distributions, and aerosol number densities differed by more than an order of magnitude for iron and stony meteoroids of the same entrance velocity. Less viscous iron melt yielded finer particles and a denser aerosol wake than the stony one. The fast iron meteoroid with the entry velocity V_∞ = 60 km s^-1 gradually burst into more than a million particles with radii ranging from 4 μm to 135 μm when the meteoroid radius was 0.3 cm. Generally, the radii of the sprayed droplets are comparable with sizes of the “grains” obtained within the theories based on the “dustball” hypothesis. However, the parameters of sprayed droplets depended essentially on a meteoroid radius, velocity and substance. A future effort may involve modelling the luminosity of the meteoroid and sprayed droplets to compare with observed meteor light curves or trail lengths.

During the spraying process the largest droplets are generated at the beginning of the ablation, while the smallest droplets are generated towards the end, as the present results show. The ablating meteoroid radius varies almost linearly with time. The influence of the meteoroid entry velocity on the number and sizes of breakaway particles is estimated.

This paper presents a quantitative application of the gradient instability theory to meteoroid fragmentation. Based on the results of this study, the next steps are investigation of the dynamic structure of the meteoric aerosol wake and meteor luminosity under the increasing rate of vapour production caused by the increased liquid surface of sprayed droplets, as well as investigating the case of an arbitrary meteoroid entry angle.

Movies

Movie of Fig. 6 Access here

Movie of Fig. 7 Access here

References

Appendix A: The cause of the gradient instability

The cause of instability of a gradient flow in a vicinity of liquid – gas interface is the inertial forces of disturbed curvilinear motion (Fig. 1c). Indeed, in its motion along curvilinear trajectory with velocity V every liquid particle undergoes centrifugal acceleration: $\begin{array}{ccc}{\mathit{a}}_{\mathrm{cf}}\mathrm{=}\frac{{\mathit{V}}^{\mathrm{2}}}{{\mathit{R}}_{\mathrm{c}}}\mathrm{\approx }\frac{{\mathit{V}}_{\mathrm{s}}^{\mathrm{2}}}{{\mathit{\delta }}_{\mathrm{m}}}\mathrm{=}\frac{\mathrm{(}\mathit{\alpha \mu }{\mathrm{)}}^{\mathrm{2}\mathit{/}\mathrm{3}}}{{\mathit{p}}^{\mathrm{2}}}\frac{{\mathit{V}}_{\mathrm{a}}^{\mathrm{2}}}{{\mathit{\delta }}_{\mathrm{a}}}\frac{{\mathit{\mu }}^{\mathrm{2}\mathit{/}\mathrm{3}}}{{\mathit{\alpha }}^{\mathrm{1}\mathit{/}\mathrm{3}}}\mathrm{=}\frac{{\mathit{\alpha }}^{\mathrm{1}\mathit{/}\mathrm{3}}{\mathit{\mu }}^{\mathrm{4}\mathit{/}\mathrm{3}}}{{\mathit{p}}^{\mathrm{2}}}\frac{{\mathit{V}}_{\mathrm{a}}^{\mathrm{2}}}{{\mathit{\delta }}_{\mathrm{a}}}\mathit{,}& & \end{array}$(A.1)where Eqs. (2) are used. Here the curvature radius R_c = AB (Fig. 1c) of the disturbed trajectory is proportional to ≈ λ_f/ 2 ≈ δ_m, as the expression Δ_f = 2πδ_m/λ_f = 1.225 shows at We_s> 20 (see Sect. 3).

The centrifugal force is normal to the air-melt interface, acting on any liquid particle within the disturbance elevation (Fig. 1c). Thus, once the melt surface is disturbed, the liquid heaps are being pulled out by the centrifugal force with huge acceleration of the order of ≈10⁴ − 10⁵ m s^-2, as Eq. (A.1) shows for meteoroids. These values are more than a thousand times greater than the terrestrial gravitation, which explains the mechanism of the ligament formation shown in Fig. 4. The same tendency has the centrifugal force acting on the moving air particles above the liquid heap, but this action at the interface is smaller because of the less air density. Indeed, the dispersal relation (Eq. (1)) of the perturbation problem contains none of the gas parameters, so, the air contributes a little force into the interface instability mechanism, its role is reduced simply to generate the shear motion in liquid by means of the viscous stresses.

The huge centrifugal force pulls out the liquid ligaments of size ≈λ/ 2, which quickly burst by the airstream into several drops of approximately the same diameter (Fig. 4), providing the mechanism of meteoroid melt spraying. The time of accelerated growth of the disturbance amplitude till the value λ/ 2, which is sufficient for the breakaway droplet formation, can be estimated as t = (λ/ 2a_cf)^0.5< 10^-4 s (when Eq. (A.1) is taken into account), assuming constant acceleration. These estimates show that the unstable disturbance develops quickly during meteoroid flight thus providing the spraying ablation mechanism.

Appendix B: The boundary layer relationships

Since an air-liquid interaction occurs in the conjugated boundary layers, a few relationships are necessary for the gradient instability analysis and its application, which have their origin in the boundary layer theory. In accordance with this theory (Schlichting 1960), the air, δ_a, and the melt, δ_m, boundary layer thicknesses are expressed via the Reynolds numbers of flows in the air, Re_a, and melt, Re_m, boundary layers: $\begin{array}{ccc}{\mathit{\delta }}_{\mathrm{a}}\mathrm{=}\mathit{G}\frac{\mathrm{2}{\mathit{R}}_{\mathrm{0}}}{{\mathrm{Re}}_{\mathrm{a}}^{\mathrm{0}\mathit{.5}}}\mathit{,}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathit{\delta }}_{\mathrm{m}}\mathrm{=}\mathit{G}\frac{\mathrm{2}{\mathit{R}}_{\mathrm{0}}}{{\mathrm{Re}}_{\mathrm{m}}^{\mathrm{0.5}}}\mathit{,}& & \end{array}$(B.1)here Re_a = ρ_a2R₀(V_a − V_s) /μ_a, Re_m = ρ_m2R₀V_s/μ_m, μ_a is the air dynamic viscosity, V_a is the air velocity, V_s is velocity at the air-melt interface (Fig. 2,b), R₀ is the initial meteoroid radius. Characteristic values of the velocity gradients in boundary layers are bounded by the requirement of equality of the viscous shear stresses at the interface y = 0: μ_mV_s/δ_m = μ_a(V_a − V_s) /δ_a. Eliminating geometry factor G, we obtain the important relations: $\begin{array}{ccc}{\mathit{\delta }}_{\mathrm{m}}\mathrm{=}{\left(\frac{\mathit{\alpha }}{{\mathit{\mu }}^{\mathrm{2}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{3}}{\mathit{\delta }}_{\mathrm{a}}\mathit{,}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathit{V}}_{\mathrm{s}}\mathrm{=}\frac{{\left(\mathit{\alpha \mu }\right)}^{\mathrm{1}\mathit{/}\mathrm{3}}}{\mathit{p}}{\mathit{V}}_{\mathrm{a}}\mathit{.}& & \end{array}$(B.2)Assuming the spherical shape of the meteoroid and the velocity potential of the air flow around a sphere in the form V_a(ϕ) = 1.5(V_∞ − w)sinϕ (Kochin et al. 1964), let us apply the boundary layer thickness distribution along the meteoroid surface in Ranger’s (1972) form: ${\mathit{\delta }}_{\mathrm{a}}\mathrm{\left(}\mathit{\varphi ,t}\mathrm{\right)}\mathrm{=}\mathrm{2.2}\mathit{R}\mathrm{\left(}\mathit{t}\mathrm{\right)}{\mathrm{Re}}_{\mathrm{a}}^{-0.5}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{\Psi }\mathrm{\left(}\mathit{\varphi }\mathrm{\right)}$; the latter is valid for a gas over a liquid sphere, here Ψ(ϕ) = ((6ϕ − 4sin2ϕ + 0.5sin4ϕ) / sin⁵ϕ)^0.5, ϕ is the polar angle of a surface area element (Figs. 1b and 3). The latter together with Eq. (B.1) yields the geometry factor for a liquid sphere: $\mathit{G}\mathrm{\left(}\mathit{\varphi ,t}\mathrm{\right)}\mathrm{=}\mathrm{1.1}\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{R}\end{array}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{\Psi }\mathrm{\left(}\mathit{\varphi }\mathrm{\right)}$, $\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{R}\end{array}\mathrm{=}\mathit{R}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathit{/}{\mathit{R}}_{\mathrm{0}}$.

Let us transform the necessary condition for the gradient instability, We_s ≥ 3.08, rewriting it for an arbitrary area element ϕ at meteoroid surface: $\begin{array}{ccc}{\mathrm{We}}_{\mathrm{s}}\mathrm{\left(}\mathit{\phi ,\tau }\mathrm{\right)}\mathrm{=}\frac{{\mathit{\rho }}_{\mathrm{m}}{\mathit{V}}_{\mathrm{s}}^{\mathrm{2}}{\mathit{\delta }}_{\mathrm{m}}}{\mathrm{\Sigma }}\mathrm{=}\frac{{\mathit{\rho }}_{\mathrm{\infty }}{\mathit{V}}_{\mathrm{\infty }}^{\mathrm{2}}\mathrm{2}{\mathit{R}}_{\mathrm{0}}}{\mathrm{\Sigma }}\frac{{\mathit{\rho }}_{\mathrm{m}}}{{\mathit{\rho }}_{\mathrm{a}}}\frac{{\mathit{\rho }}_{\mathrm{a}}}{{\mathit{\rho }}_{\mathrm{\infty }}}\frac{{\mathit{V}}_{\mathrm{s}}^{\mathrm{2}}}{{\mathit{V}}_{\mathrm{a}}^{\mathrm{2}}}\frac{{\mathit{V}}_{\mathrm{a}}^{\mathrm{2}}}{{\mathit{V}}_{\mathrm{\infty }}^{\mathrm{2}}}\frac{{\mathit{\delta }}_{\mathrm{m}}}{{\mathit{\delta }}_{\mathrm{a}}}\frac{{\mathit{\delta }}_{\mathrm{a}}}{\mathrm{2}{\mathit{R}}_{\mathrm{0}}}\mathrm{·}& & \end{array}$(B.3)The right-hand side can be further transformed to: $\begin{array}{ccc}& & \mathrm{=}{\mathrm{We}}_{\mathrm{\infty }}\frac{\mathrm{1}}{\mathit{\alpha }}\frac{{\mathit{\rho }}_{\mathrm{a}}}{{\mathit{\rho }}_{\mathrm{\infty }}}\frac{{\mathit{V}}_{\mathrm{a}}^{\mathrm{2}}}{{\mathit{V}}_{\mathrm{\infty }}^{\mathrm{2}}}\frac{\mathrm{(}\mathit{\alpha \mu }{\mathrm{)}}^{\mathrm{2}\mathit{/}\mathrm{3}}}{{\mathit{p}}^{\mathrm{2}}}\frac{{\mathit{\alpha }}^{\mathrm{1}\mathit{/}\mathrm{3}}}{{\mathit{\mu }}^{\mathrm{2}\mathit{/}\mathrm{3}}}\mathit{G}\frac{\begin{array}{c}\mathrm{˜}\\ \mathit{R}\end{array}}{\mathrm{R}{\mathrm{e}}_{\mathrm{\infty }}^{\mathrm{0.5}}}\frac{\mathrm{R}{\mathrm{e}}_{\mathrm{\infty }}^{\mathrm{0.5}}}{\mathrm{R}{\mathrm{e}}_{\mathrm{a}}^{\mathrm{0.5}}}\\ & & \mathrm{=}\frac{{\mathrm{We}}_{\mathrm{\infty }}}{\mathrm{R}{\mathrm{e}}_{\mathrm{\infty }}^{\mathrm{0.5}}}\frac{{\mathit{\rho }}_{\mathrm{a}}}{{\mathit{\rho }}_{\mathrm{\infty }}}\frac{{\mathit{V}}_{\mathrm{a}}^{\mathrm{2}}}{{\mathit{V}}_{\mathrm{\infty }}^{\mathrm{2}}}\frac{\mathit{G}}{{\mathit{p}}^{\mathrm{2}}}\frac{\begin{array}{c}\mathrm{˜}\\ \mathit{R}\end{array}}{\begin{array}{c}\mathrm{˜}\\ {\mathit{R}}^{\mathrm{0.5}}\end{array}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{W}{\mathrm{)}}^{\mathrm{0.5}}}\mathrm{·}\end{array}$(B.4)When ${\mathit{\rho }}_{\mathrm{a}}\mathrm{=}{\mathit{\rho }}_{\mathrm{\infty }}\mathit{,}{\mathit{V}}_{\mathrm{a}}^{\mathrm{2}}\mathrm{=}{\mathit{V}}_{\mathrm{\infty }}^{\mathrm{2}}$, and in the virtue of relationships $\begin{array}{ccc}{\mathit{\delta }}_{\mathrm{m}}\mathrm{=}{\left(\frac{\mathit{\alpha }}{{\mathit{\mu }}^{\mathrm{2}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{3}}{\mathit{\delta }}_{\mathrm{a}}\mathrm{;}\hspace{0.17em}\hspace{0.17em}{\mathit{V}}_{\mathrm{s}}\mathrm{=}\frac{{\left(\mathit{\alpha \mu }\right)}^{\mathrm{1}\mathit{/}\mathrm{3}}}{\mathit{p}}{\mathit{V}}_{\mathrm{a}}\mathrm{;}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathit{\delta }}_{\mathrm{a}}\mathrm{=}\mathit{G}\frac{\mathrm{2}\mathit{R}\mathrm{\left(}\mathit{t}\mathrm{\right)}}{{\mathrm{R}{\mathrm{e}}_{\mathrm{a}}}^{\mathrm{0.5}}\mathrm{\left(}\mathit{t}\mathrm{\right)}}\mathrm{;}& & \end{array}$V_a = 1.5V_∞(1 − w(t))sin(ϕ), we will have: $\begin{array}{ccc}\mathrm{R}{\mathrm{e}}_{\mathrm{a}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\mathrm{R}{\mathrm{e}}_{\mathrm{\infty }}\frac{\mathit{R}\mathrm{\left(}\mathit{t}\mathrm{\right)}}{{\mathit{R}}_{\mathrm{0}}}\frac{{\mathit{V}}_{\mathrm{\infty }}\mathrm{-}\mathit{w}\mathrm{\left(}\mathit{t}\mathrm{\right)}}{{\mathit{V}}_{\mathrm{\infty }}}\mathrm{=}\mathrm{R}{\mathrm{e}}_{\mathrm{\infty }}\begin{array}{c}\mathrm{˜}\\ \mathit{R}\end{array}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{W}\mathrm{(}\mathit{\tau }\mathrm{\left)}\mathrm{\right)}\mathit{.}& & \end{array}$Eventually, we obtain Eq. (3): ${\mathrm{We}}_{\mathrm{s}}\mathrm{\left(}\mathit{\varphi ,\tau }\mathrm{\right)}\mathrm{=}\frac{\mathrm{2.475}}{{\mathit{p}}^{\mathrm{2}}}\sqrt{\begin{array}{c}\mathrm{˜}\\ \mathit{R}\end{array}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}{\left(\mathrm{1}\mathrm{-}\mathit{W}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\right)}^{\mathrm{3}}}{\sin }^{\mathrm{2}}\mathit{\varphi }\mathrm{\Psi }\mathrm{\left(}\mathit{\varphi }\mathrm{\right)}\hspace{0.17em}\mathrm{GI}\mathrm{\ge }\mathrm{3.08.}$(B.5)Thus, the local condition for the gradient instability proceeding, We_s> We_s.cr. = 3.08, can be rewritten in the whole problem velocity and length scales, V_∞ and R₀: $\begin{array}{ccc}\mathrm{GI}\mathrm{=}{\mathrm{We}}_{\mathrm{\infty }}{\mathrm{Re}}_{\mathrm{\infty }}^{-0.5}\hspace{0.17em}\hspace{0.17em}\mathit{>}{\mathrm{GI}}_{\mathrm{cr}}\hspace{0.17em}\mathrm{=}\mathrm{0}\mathit{.4}\mathit{,}& & \end{array}$(B.6)which exhibits the necessary condition for the melt surface instability for the whole meteoroid flight. Expression (B.6) is obtained from the inequality in (B.5) at τ = 0 since $\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{R}\end{array}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}$, 1 − W(τ) are the decreasing functions, so that initial moment is preferable for instability performance. Value ϕ = π/ 2 was also used since δ_m,V_s, We_s have the uppermost values at a liquid body periphery, ϕ = π/ 2.

Appendix C: The meaning of the surface Weber number criterion

The two Weber numbers, surface We_s and meteoroid We_∞, have the same general sense: the ratio of the hydrodynamic forces to the surface tension forces. But the surface Weber number We_s specifically relates hydrodynamic pressure ${\mathit{\rho }}_{\mathrm{m}}{\mathit{V}}_{\mathrm{s}}^{\mathrm{2}}\mathit{/}\mathrm{2}$ of the melt motion in the liquid boundary layer to the surface tension forces δ_m/ Σ of the disturbed interface, whose curvature radius R_c = AB (Fig. 1c) is proportional to δ_m; the latter follows from the relation Δ_f = 2πδ_m/λ_f = 1.225 at We_s> 20. The former force creates the destabilizing effect, defining the centrifugal force scale, see Eq. (A.1); the latter force defines a stabilizing effect. Thus, the surface Weber number reflects locally, at a disturbance level, the ratio of destabilizing forces to the stabilizing ones, governing the disturbance development. At We_s> 3.08 the inertial forces prevail over the surface tension and the disturbance becomes unstable, growing in amplitude. The role of viscous forces is to form the gradient flow of characteristic depth δ_m in a liquid melt.

The existence of the critical value We_s.cr. has a simple hydrodynamic sense. In order to reveal the meaning of We_s.cr for a meteoroid flight, let us go over to the air flow parameters in the inequality $\begin{array}{ccc}\frac{{\mathit{\rho }}_{\mathrm{m}}{\mathit{V}}_{\mathrm{s}}^{\mathrm{2}}{\mathit{\delta }}_{\mathrm{m}}}{\mathrm{\Sigma }}\mathrm{=}\hspace{0.17em}\frac{{\mathit{\rho }}_{\mathrm{a}}{\mathit{V}}_{\mathrm{a}}^{\mathrm{2}}{\mathit{\delta }}_{\mathrm{a}}}{\mathrm{\Sigma }}\frac{{\mathit{\rho }}_{\mathrm{m}}}{{\mathit{\rho }}_{\mathrm{a}}}\frac{{\mathit{V}}_{\mathrm{s}}^{\mathrm{2}}}{{\mathit{V}}_{\mathrm{a}}^{\mathrm{2}}}\frac{{\mathit{\delta }}_{\mathrm{m}}}{{\mathit{\delta }}_{\mathrm{a}}}\mathit{>}\mathrm{3.08}\mathit{,}& & \end{array}$(C.1)or, using Eqs. (2) (or equally (B.2)): $\begin{array}{ccc}\frac{{\mathit{\rho }}_{\mathrm{a}}{\mathit{V}}_{\mathrm{a}}^{\mathrm{2}}{\mathit{\delta }}_{\mathrm{a}}}{\mathrm{\Sigma }}\frac{\mathrm{1}}{\mathit{\alpha }}\frac{\mathrm{(}\mathit{\alpha \mu }{\mathrm{)}}^{\mathrm{2}\mathit{/}\mathrm{3}}}{{\mathit{p}}^{\mathrm{2}}}\frac{{\mathit{\alpha }}^{\mathrm{1}\mathit{/}\mathrm{3}}}{{\mathit{\mu }}^{\mathrm{2}\mathit{/}\mathrm{3}}}\mathit{>}\mathrm{3.08.}& & \end{array}$(C.2)The inequality (C.2) can be resolved with respect to the dynamic pressure P of the airstream past meteoroid: $\begin{array}{ccc}\mathit{P}\mathrm{=}\frac{{\mathit{\rho }}_{\mathrm{a}}{\mathit{V}}_{\mathrm{a}}^{\mathrm{2}}}{\mathrm{2}}\mathit{>}\frac{\mathrm{3.08}}{\mathrm{2}}\mathrm{·}\frac{\mathit{\alpha }{\mathit{p}}^{\mathrm{2}}}{\mathrm{(}\mathit{\alpha \mu }{\mathrm{)}}^{\mathrm{2}\mathit{/}\mathrm{3}}}\frac{{\mathit{\mu }}^{\mathrm{2}\mathit{/}\mathrm{3}}}{{\mathit{\alpha }}^{\mathrm{1}\mathit{/}\mathrm{3}}}\frac{\mathrm{\Sigma }}{{\mathit{\delta }}_{\mathrm{a}}}\mathrm{=}\mathrm{1.54}{\mathit{p}}^{\mathrm{2}}\frac{\mathrm{\Sigma }{\mathrm{Re}}_{\mathrm{a}}^{\mathrm{0}\mathit{.5}}}{\mathrm{2}{\mathit{R}}_{\mathrm{0}}\mathit{G}}\mathit{,}& & \end{array}$(C.3)with the use of Eqs. (B.1). With the value G(π/ 2,0) = 3.38, which follows from the Ranger expression for δ_a (see Appendix B), we consecutively obtain the values for the three variants of meteoroids: P^I> 3.0 kPa, P^II> 2.1 kPa, P^III> 0.9 kPa, the notation in Sect. 3 regarding the value We_s.cr. = 4.62 being taken into account.

In some publications (Borovička et al. 2007) it has been noted that meteoroid fragmentation can not be explained by the direct action of the dynamic pressure of airstream. The latter has an order of (5 − 20) kPa at heights of about 90–100 km, which is too small to prevail over the solid meteoroid strength. But the critical values of the dynamic pressure P, obtained above, which correspond to the value We_s.cr. = 4.62, show that the dynamic pressure of airstream at heights up to 100 km, which is too low with respect to a solid meteoroid body strength, is sufficient, however, to overcome surface tension forces and to pull out liquid particles from the meteoroid melt by means of unstable disturbances, providing thus the quasi-continuous spraying mechanism of meteoroid ablation. The gradient instability pulls out ligaments from the liquid (Fig. 4) by means of centrifugal forces imposed at moving liquid layer, when the unstable disturbances appear on a liquid surface (Fig. 1c).

Appendix D: The differential equation for the sprayed droplets’ number

The rate of the tori production (see Sect. 4) can be expressed as the ratio of the number n_tor of the tori that were broken away by the time period t_per, to this time period: $\begin{array}{ccc}\mathrm{\Delta }\mathit{ṅ}\mathrm{tor}\mathrm{=}\frac{{\mathit{n}}_{\mathrm{tor}}}{{\mathit{t}}_{\mathrm{per}}}\mathrm{=}\frac{\mathit{R}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{Im}\mathrm{\left(}{\mathrm{\Omega }}_{\mathrm{f}}\mathrm{\right)}{\mathit{V}}_{\mathrm{s}}{\mathrm{\Delta }}_{\mathrm{f}}}{\mathrm{2}\mathit{\pi }{\mathit{\delta }}_{\mathrm{m}}{\mathit{k}}_{\mathrm{t}}{\mathit{\delta }}_{\mathrm{m}}}\mathrm{\Delta }\mathit{\varphi }& & \end{array}$$\begin{array}{ccc} \mathrm{=}\frac{\mathrm{Im}\mathrm{\left(}{\mathrm{\Omega }}_{\mathrm{f}}\mathrm{\right)}{\mathrm{\Delta }}_{\mathrm{f}}\mathrm{\left(}\mathit{\alpha \mu }{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{3}}\mathit{R}\mathrm{\right(}\mathit{t}\mathrm{\left)}\mathrm{1.5}{\mathit{V}}_{\mathrm{\infty }}\mathrm{\right(}\mathrm{1}\mathrm{-}\mathit{w}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{)}\sin \mathit{\varphi }}{{\mathit{k}}_{\mathrm{t}}\mathit{p}\mathrm{2}\mathit{\pi }\mathrm{(}\mathit{\alpha }\mathit{/}{\mathit{\mu }}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{3}}{\mathit{G}}^{\mathrm{2}}\mathrm{4}\mathit{R}\mathrm{(}\mathit{t}\mathrm{)}\mathit{/}\mathrm{R}{\mathrm{e}}_{\mathrm{a}}\mathrm{(}\mathit{t}\mathrm{)}}\mathrm{\Delta }\mathit{\varphi ,}& & \end{array}$(D.1)with the use of the expressions $\begin{array}{ccc}& {\mathit{\lambda }}_{\mathrm{f}}\mathrm{=}\mathrm{2}\mathit{\pi }{\mathit{\delta }}_{\mathrm{m}}\mathit{/}{\mathrm{\Delta }}_{\mathrm{f}}\mathrm{;}\hspace{0.17em}\hspace{0.17em}{\mathit{\delta }}_{\mathrm{m}}\mathrm{=}\mathrm{(}\mathit{\alpha }\mathit{/}{\mathit{\mu }}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{3}}{\mathit{\delta }}_{\mathrm{a}}\mathrm{;}\hspace{0.17em}\hspace{0.17em}{\mathit{\delta }}_{\mathrm{a}}\mathrm{=}\mathit{G}\mathrm{2}\mathit{R}\mathrm{(}\mathit{t}\mathrm{\left)}{\mathrm{R}{\mathrm{e}}_{\mathrm{a}}}^{-0.5}\mathrm{\right(}\mathit{t}\mathrm{)}\mathrm{;}& \\ & {\mathit{V}}_{\mathrm{a}}\mathrm{=}\mathrm{1.5}\mathrm{(}{\mathit{V}}_{\mathrm{\infty }}\mathrm{-}\mathit{w}\mathrm{)}\sin \mathrm{\left(}\mathit{\varphi }\mathrm{\right)}\mathrm{;}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathit{t}}_{\mathrm{per}}\mathrm{=}{\mathit{k}}_{\mathrm{t}}{\mathrm{Im}}^{-1}\mathrm{\left(}{\mathrm{\Omega }}_{\mathrm{f}}\mathrm{\right)}{\mathit{\delta }}_{\mathrm{m}}\mathit{/}{\mathrm{V}}_{\mathrm{s}}\mathrm{;}& \\ & \mathrm{R}{\mathrm{e}}_{\mathrm{a}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\mathrm{R}{\mathrm{e}}_{\mathrm{\infty }}\frac{\mathit{R}\mathrm{\left(}\mathit{t}\mathrm{\right)}}{{\mathit{R}}_{\mathrm{0}}}\frac{{\mathit{V}}_{\mathrm{\infty }}\mathrm{-}\mathit{w}\mathrm{\left(}\mathit{t}\mathrm{\right)}}{{\mathit{V}}_{\mathrm{\infty }}}\mathrm{=}\mathrm{R}{\mathrm{e}}_{\mathrm{\infty }}\begin{array}{c}\mathrm{˜}\\ \mathit{R}\end{array}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{W}\mathrm{(}\mathit{\tau }\mathrm{\left)}\mathrm{\right)}\mathit{.}& \end{array}$The single torus volume is υ_tor = 2π²(k_rλ_f)²R(t)sinϕ, so that the melt volume gone with the tori by the time period t_per is $\mathrm{\Delta }{\mathit{\upsilon }}_{\mathrm{tor}}\mathrm{=}\mathrm{\Delta }\mathit{ṅ}{}_{\mathrm{tor}}\mathit{\upsilon }_{\mathrm{tor}}$. Assuming that the liquid tori are quickly breaking up into droplets of radius r in the airstream and dividing this volume by the single droplet volume, υ_dr = 4 / 3π(k_rλ_f)³, we obtain the number of droplets broken away by the time t_per: $\begin{array}{ccc}\mathrm{\Delta }\mathit{ṅ}\mathrm{dr}& \mathrm{=}& \frac{\mathrm{\Delta }\mathit{ṅ}\mathrm{tor}{\mathit{\upsilon }}_{\mathrm{tor}}}{{\mathit{\upsilon }}_{\mathrm{dr}}}\\ & \mathrm{=}& \frac{\mathrm{Im}\mathrm{\left(}{\mathrm{\Omega }}_{\mathrm{f}}\mathrm{\right)}{\mathrm{\Delta }}_{\mathrm{f}}\mathrm{(}\mathit{\alpha \mu }{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{3}}\mathit{R}\mathrm{·}\mathrm{1.5}\mathrm{(}{\mathit{V}}_{\mathrm{\infty }}\mathrm{-}\mathit{w}\mathrm{\left)}\sin \mathit{\varphi }\mathrm{2}{\mathit{\pi }}^{\mathrm{2}}\mathrm{\right(}{\mathit{k}}_{\mathrm{r}}{\mathit{\lambda }}_{\mathrm{f}}{\mathrm{)}}^{\mathrm{2}}\mathit{R}\mathrm{·}\sin \mathit{\varphi }}{{\mathit{k}}_{\mathrm{t}}\mathit{p}\mathrm{2}\mathit{\pi }\mathrm{(}\mathit{\alpha }\mathit{/}{\mathit{\mu }}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{3}}{\mathit{G}}^{\mathrm{2}}\mathrm{4}\mathit{R}\mathit{/}\mathrm{R}{\mathrm{e}}_{\mathrm{a}}\mathrm{·}\mathrm{4}\mathit{/}\mathrm{3}\mathit{\pi }\mathrm{(}{\mathit{k}}_{\mathrm{r}}{\mathit{\lambda }}_{\mathrm{f}}{\mathrm{)}}^{\mathrm{3}}}\mathrm{\Delta }\mathit{\varphi }\mathit{.}\\ & & \end{array}$(D.2)Using the above expressions for λ_f, δ_m, δ_a, Re_a, and G = 1.1Ψ(ϕ), we can simplify this expression to the form: $\begin{array}{ccc}\mathrm{\Delta }\mathit{n}\mathrm{\left(}\mathit{\varphi ,\tau }\mathrm{\right)}\mathrm{=}{\mathit{B}}_{\mathrm{2}}\sqrt{\begin{array}{c}\mathrm{˜}\\ \mathit{R}\end{array}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}{\left(\mathrm{1}\mathrm{-}\mathit{W}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\right)}^{\mathrm{5}}}\frac{{\sin }^{\mathrm{2}}\mathit{\varphi }}{{\mathrm{\Psi }}^{\mathrm{3}}\mathrm{\left(}\mathit{\varphi }\mathrm{\right)}}\mathrm{\Delta }\mathit{\varphi }\mathrm{\Delta }\mathit{\tau }\mathit{.}& & \end{array}$(D.3)Here ${\mathit{B}}_{\mathrm{2}}\mathrm{=}\frac{\mathrm{0.21}{\mathrm{\Delta }}_{\mathrm{f}}^{\mathrm{2}}\mathrm{Im}\mathrm{\left(}{\mathrm{\Omega }}_{\mathrm{f}}\mathrm{\right)}\sqrt{\mathrm{2}{\mathrm{Re}}_{\mathrm{\infty }}^{\mathrm{3}}}}{\mathit{\pi }\hspace{0.17em}{\mathit{k}}_{\mathrm{r}}{\mathit{k}}_{\mathrm{t}}\mathrm{(}\mathrm{1}\mathrm{+}\mathrm{(}\mathit{\alpha \mu }{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{3}}\mathrm{)}}\frac{{\mathit{\mu }}^{\mathrm{7}\mathit{/}\mathrm{3}}}{{\mathit{\alpha }}^{\mathrm{7}\mathit{/}\mathrm{6}}}$ has the sense of the scaling parameter for the number of breakaway particles (Girin 2011a).

Equation (D.3) is the differential Eq. (7) for the number Δn of the newborn particles, which are released during the period Δτ = t_per/t_ch from the area element Δl = R(t)Δϕ.

In accord with Eq. (5) and the expressions given above, the newborn particle radius can be written as follows: $\begin{array}{ccc}\mathit{r}\mathrm{\left(}\mathit{\varphi ,\tau }\mathrm{\right)}\mathrm{=}{\mathit{k}}_{\mathrm{r}}\frac{\mathrm{2}\mathit{\pi }{\mathit{\delta }}_{\mathrm{m}}}{{\mathrm{\Delta }}_{\mathrm{f}}}\mathrm{=}{\mathit{k}}_{\mathrm{r}}\frac{\mathrm{2}\mathit{\pi }\mathrm{·}\mathrm{2}{\mathit{R}}_{\mathrm{0}}\mathrm{·}\mathrm{1.1}\mathrm{\Psi }\mathrm{\left(}\mathit{\varphi }\mathrm{\right)}}{{\mathrm{\Delta }}_{\mathrm{f}}{\mathrm{R}{\mathrm{e}}_{\mathrm{\infty }}}^{\mathrm{0.5}}}\frac{\begin{array}{c}\mathrm{˜}\\ \mathit{R}\end{array}\mathrm{(}\mathit{\tau }{\mathrm{)}}^{\mathrm{0.5}}}{\mathrm{[}\mathrm{1}\mathrm{-}\mathit{W}\mathrm{(}\mathit{\tau }\mathrm{)}{\mathrm{]}}^{\mathrm{0.5}}}\frac{{\mathit{\alpha }}^{\mathrm{1}\mathit{/}\mathrm{3}}}{{\mathit{\mu }}^{\mathrm{2}\mathit{/}\mathrm{3}}}\mathit{.}& & \end{array}$(D.4)This expression can be rewritten in the form of Eq. (8): $\begin{array}{ccc}\mathit{r}\mathrm{\left(}\mathit{\varphi ,\tau }\mathrm{\right)}\mathrm{=}{\mathit{R}}_{\mathrm{0}}{\mathit{B}}_{\mathrm{1}}\hspace{0.17em}\sqrt{\frac{\begin{array}{c}\mathrm{˜}\\ \mathit{R}\end{array}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}}{\mathrm{1}\mathrm{-}\mathit{W}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}}}\mathrm{\Psi }\mathrm{\left(}\mathit{\varphi }\mathrm{\right)}\mathit{.}& & \end{array}$(D.5)Here ${\mathit{B}}_{\mathrm{1}}\mathrm{=}\frac{\mathrm{4.4}\mathit{\pi }\hspace{0.17em}{\mathit{k}}_{\mathrm{r}}}{{\mathrm{\Delta }}_{\mathrm{f}}\sqrt{\mathrm{2}{\mathrm{Re}}_{\mathrm{\infty }}}}\frac{{\mathit{\alpha }}^{\mathrm{1}\mathit{/}\mathrm{3}}}{{\mathit{\mu }}^{\mathrm{2}\mathit{/}\mathrm{3}}}$ has  the  sense  of  the  scaling  parameter for the breakaway particle sizes (Girin 2011a). For large values of We_s, when We_s> 20, B₁ and B₂ are constants since Δ_f and Ω_f are constants (Fig. 5).

Appendix E: Nomenclature

b	– molten layer thickness;
B 1 , B 2	– scaling parameters;
E	– disturbance amplitude;
G	– boundary layer geometry factor;
$\mathrm{GI}\mathrm{=}{\mathrm{We}}_{\mathrm{\infty }}{\mathrm{Re}}_{\mathrm{\infty }}^{-0.5}$	– gradient instability criterion;
Im-1(ω)	– growth rate factor;
k r , k t	– coefficients in Eqs. (5) and (6);
l	– meteoroid path;
m	– meteoroid mass;
M = m/m0
Ma	– Mach number;
N	– total number of particles;
p = 1 + (αμ)1 / 3
P	– dynamic pressure;
R	– current meteoroid radius;
$\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{R}\end{array}\mathrm{=}\mathit{R}\mathit{/}{\mathit{R}}_{\mathrm{0}}$
r	– breakaway droplet radius;
$\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{r}\end{array}\mathrm{=}\mathit{r}\mathit{/}{\mathit{R}}_{\mathrm{0}}$
Re∞ = ρ∞2R0V∞/μ∞	– Reynolds number for meteoroid;
S	– meteoroid surface area;
t	– elapsed time;
tch = 2R0/α0.5V∞	– characteristic timescale;
V	– velocity;
w	– meteoroid velocity deficit due to
	deceleration;
W = w/V∞
${\mathrm{We}}_{\mathrm{\infty }}\mathrm{=}{\mathit{\rho }}_{\mathrm{\infty }}{\mathit{V}}_{\mathrm{\infty }}^{\mathrm{2}}\mathrm{2}{\mathit{R}}_{\mathrm{0}}\mathit{/}\mathrm{\Sigma }$	– Weber number for meteoroid;
${\mathrm{We}}_{\mathrm{s}}\mathrm{=}{\mathit{\rho }}_{\mathrm{m}}{\mathit{V}}_{\mathrm{s}}^{\mathrm{2}}{\mathit{\delta }}_{\mathrm{m}}\mathit{/}\mathrm{\Sigma }$	– surface Weber number;
x,y	– coordinates within boundary layer;
z	– meteor altitude.

All Tables

All Figures

	Fig. 1 Scheme of the flow in case of outstripping ablation (a) and in case of outstripping melting (b). Scheme of flow in disturbed boundary layer (c): solid thin lines indicate melt particles’ streamlines; thick dashed lines indicate undisturbed boundary layer limits. AB is the curvature radius of disturbed streamlines; Fcf is the centrifugal force acting at moving melt particle.
In the text

	Fig. 2 Linear fit of velocity profile in conjugated boundary layers; case νm>νa(a) and νa>νm(b); δa, δm – air and melt boundary layer thicknesses.
In the text

	Fig. 3 Flow scheme around meteoroid body. Stable and unstable parts of molten meteoroid surface are shown.
In the text

	Fig. 4 Particle release process at ith area element of meteoroid surface.
In the text

	Fig. 5 Basic dependencies Δf(Wes) (solid line) and Im [ Ωf(Wes) ] (dashed line) of gradient instability.
In the text

	Fig. 6 Distributions of breakaway droplets by sizes: Δn(r)(a) and ΔM(r)(b). Variant I: fast iron meteoroid. Blue indicates the final distribution at τ = 2.72; brown, red, light green, and crimson indicate, respectively, the intermediate distributions at τ = 0.44, 0.89, 1.39, and  1.80; tch = 2.16 × 10-3 s. See Movie 1.
In the text

	Fig. 7 Distributions of breakaway droplets by sizes: Δn(r)(a) and ΔM(r)(b). Variant II: slow iron meteoroid. Blue indicates the final distribution at τ = 3.11; brown, red, light green, and crimson indicate, respectively, the intermediate distributions at τ = 0.30, 0.93, 1.88, and 2.48; tch = 5.15 × 10-3 s. See Movie 2.
In the text

	Fig. 8 Distributions of breakaway droplets by sizes: Δn(r)(a) and ΔM(r)(b). Variant III: stony meteoroid. Blue indicates the final distribution at τ = 105; purple, orange, light-blue, green, brown, red, light green, and crimson indicate, respectively, the intermediate ones at τ = 2.05, 3.88, 5.94, 9.96, 16.8, 29.9, 48.7, and 73.0;tch = 1.45 × 10-3 s.
In the text