A chemical model of meteoric ablation

Most of the extraterrestrial dust entering the Earth’s atmosphere ablates to produce metal vapours, which have significant effects on the aeronomy of the upper mesosphere and lower thermosphere. A new Chemical Ablation Model (CAMOD) is described which treats the physics and chemistry of ablation, by including the following processes: sputtering by inelastic collisions with air molecules before the meteoroid melts; evaporation of atoms and oxides from the molten particle; diffusion-controlled migration of the volatile constituents (Na and K) through the molten particle; and impact ionization of the ablated fragments by hyperthermal collisions with air molecules. Evaporation is based on thermodynamic equilibrium in the molten meteoroid (treated as a melt of metal oxides), and between the particle and surrounding vapour phase. The loss rate of each element is then determined assuming Langmuir evaporation. CAMOD successfully predicts the meteor head echo appearance heights, observed from incoherent scatter radars, over a wide range of meteoroid velocities. The model also confirms that differential ablation explains common-volume lidar observations of K, Ca and Ca + in fresh meteor trails. CAMOD is then used to calculate the injection rates into the atmosphere of a variety of elements as a function of altitude, integrated over the meteoroid mass and velocity distributions. The most abundant elements (Fe, Mg and Si) have peak injection rates around 85 km, with Na and K about 8 km higher. The more refractory element Ca ablates around 82 km with a Na:Ca ratio of 4:1, which does therefore not explain the depletion of atomic Ca to Na, by more than 2 orders of magnitude, in the upper mesosphere. Diffusion Correspondence to: J. M. C. Plane (j.m.c.plane@leeds.ac.uk) of the most volatile elements (Na and K) does not appear to be rate-limiting except in the fastest meteoroids. Nonthermal sputtering causes ∼35% mass loss from the fastest (∼60–70 km s−1) and smallest (10 −17–10−13 g) meteoroids, but makes a minor contribution to the overall ablation rate.


Introduction
There are two principal sources of meteoroids in the Earth's atmosphere (Ceplecha et al., 1998;Williams, 2002).First are the dust trails produced by sublimating comets as they orbit the sun, which are the origin of meteor showers such as the Perseids and Leonids.Second are fragments from the asteroid belt beyond Mars, and dust particles from longdecayed cometary trails.These give rise to the continuous input of sporadic meteoroids which provides a much greater mass flux than the showers (Ceplecha et al., 1998;Williams, 2002).Because of their very high entry velocities, meteoroids undergo rapid frictional heating by collision with air molecules, and their constituent minerals subsequently vaporize.This provides the dominant source of various metals and silicon in the upper atmosphere, which are manifest as layers of neutral metal atoms (Na, Fe, Ca etc.) between about 80 and 105 km, and as sporadic E layers between 90 and 140 km (Plane, 2003).Below 85 km the metals form compounds such as oxides and carbonates which condense to form meteoric smoke particles (e.g.Hunten et al., 1980;Kalashnikova et al., 2000;Saunders and Plane, 2006).These nm-sized particles most probably provide ice nuclei for the formation of noctilucent clouds in the summer high latitude mesosphere (Megner et al., 2006), and polar stratospheric clouds in the wintertime polar stratosphere (Curtius et al., 2005).
T. Vondrak et al.: Meteoric ablation The input flux of meteoroids into the atmosphere is rather uncertain, because no single technique can observe particles over the mass range from about 10 −12 to 1 g, which make up the bulk of the incoming material (Ceplecha et al., 1998).In addition, meteoroids which originate within the solar system (the vast majority) have a large range of entry velocities, from 11.5 km s −1 for a particle in the same orbit as the Earth, to 72.5 km s −1 for a particle in a retrograde orbit (Baggaley, 2002).Until the 1990s, a commonly assumed figure for the daily mass input was 44 t d −1 , averaged over the entire Earth (Hughes, 1978).This estimate was made by extrapolating between visual meteors (mass >10 mg) and satellite impact data (mass <1 µg), and ignoring indications of a somewhat lower flux from measurements made by conventional meteor radars (these instruments measure specular reflections from the ion trail created in the atmosphere by the ablating meteoroid) (Hughes, 1997).More recently, the analysis of small particle impact craters on the Long Duration Exposure Facility (LDEF), an orbital impact detector placed on a spacecraft for several years, has yielded an estimate of 110±55 t d −1 (Love and Brownlee, 1993;McBride et al., 1999).The median mass of the incoming dust particles was deduced to be ∼10 µg with a mean entry velocity of about 18 km s −1 .However, it is important to note that the LDEF experiment measured crater size, which was treated as a proxy for particle kinetic energy.Hence, the particle velocity distribution had to be assumed in order to produce the mass distribution.
In the last few years, high-powered large aperture (HPLA) radars have reported direct observations of the meteor head echo (i.e. the ball of plasma around the ablating particle as it descends through the atmosphere).This enables measurements of the direction of the line-of-sight velocity, deceleration and meteoroid mass to be made (e.g.Pellinen-Wannberg and Wannberg, 1994;Close et al., 2000;Janches et al., 2000).The observed mass distribution is shifted to smaller mass ranges (median mass ≈1 µg) compared to the LDEF results, and the total incoming mass is about 1 order of magnitude lower than the LDEF estimate.The mean entry velocity is around 40-50 km s −1 (Mathews et al., 2001), significantly higher than distributions measured by meteor radars which detect the reflection from a semi-stationary ionized trail left behind the meteoroid path (see, for example, Galligan and Baggaley (2004), Fig. 27a).The difference in the velocity distributions measured by the two techniques is due to the echo ceiling effect inherently present in the meteor radar observations.That is, these radars do not efficiently detect meteors which occur at higher altitudes (>100 km), because of the rapid diffusion of the ionized trails (Chau et al., 2007).Since faster meteors generally occur at higher altitudes, distributions measured by meteor radars are biased towards the lower speeds.In fact, Janches et al. (2008) showed that HPLA radars observe the same population of meteors as observed by meteor radars, and in addition detect a population of faster meteors that ablate at altitudes where specular trails are not efficiently detected.However, the magnitude of the head echo still depends on the meteoroid mass and velocity, and each HPLA radar is sensitive to a particular mass range (Janches et al., 2008).This implies that the velocity distribution of the smallest particles measured by an HPLA radar will be biased towards faster speeds.In particular, for the case of the Arecibo HPLA radar this effect is determined by the ablation limit (Fentzke and Janches, 2008), where small and slow particles will not have sufficient kinetic energy to ablate, and hence will not produce sufficient electrons to be detected.
The mass flux has also been estimated in other ways.Single-particle analyses of stratospheric aerosol have shown that half of the particles in the lower stratosphere contain 0.5 to 1.0 weight percent meteoric iron by mass, requiring a total extraterrestrial flux of between 20 and 100 t d −1 (Cziczo et al., 2001).Measurements of the accumulation of Ir and Pt in polar ice cores have been used to deduce a mass flux of 40±16 t d −1 (Gabrielli et al., 2004); measurements of supermagnetism from cosmic Fe in ice cores indicate a flux of 35±10 t d −1 (Lanci et al., 2007).Finally, a model of the mesospheric Na layer requires a flux of 10-30 t d −1 to match lidar observations of the layer (Plane, 2004).A consensus around 20-40 t d −1 has therefore emerged in the last 5 years, but this remains highly uncertain and is in any case a crude average that does not take account of seasonal and latitudinal variations (Janches et al., 2006;Fentzke and Janches, 2008).
Because of their very high entry velocities, meteoroids undergo rapid frictional heating by collision with air molecules, and their constituent minerals subsequently vaporize.The physics of this process has been treated in detail by several investigators.The frictional heating is balanced by radiative losses and by the absorption of heat energy through temperature increases, melting, phase transitions, and vaporization.In order to calculate these terms, parameters such as the meteoroid shape, density, and composition are needed.The question of composition has been discussed in detail recently (Rietmeijer, 2000(Rietmeijer, , 2002)).There is some uncertainty because of the great variability in composition of the different types of meteorites (Rietmeijer, 2000).Furthermore, it may be that the composition of the meteoroids that ablate in the upper atmosphere is different from that of the meteorites that survive transit through the atmosphere.Nevertheless, the current working assumption is that most of the extra-terrestrial material has the composition of ordinary chondrites.
As a starting point, most of the models developed for estimating the input of meteoric metals into the atmosphere assume that the relative metallic abundances in the ablated vapour are given by their meteoritic abundances (Plane, 1991;McNeil et al., 1995).The major metallic constituents by elemental atomic abundance are then: Mg 14.4%, Si 13.6%, Fe 12.1%, Al 1.2%, Ca 0.82% and Na 0.80% (Mason, 1971;Fegley and Cameron, 1987;Lodders and Fegley, 1988;Sears and Dodd, 1988).A more realistic approach involving fractionation has been proposed recently, with the relatively volatile elements such as Na and K evaporating first and refractory elements such as Ca evaporating last (McNeil et al., 1998).This differential ablation model, which is based on fractionation models used by planetary scientists (Fegley and Cameron, 1987), was able to explain the large depletion of Ca to Na observed in meteor trails (von Zahn et al., 1999).However, it should be borne in mind that planetary fractionation models are based on equilibrium thermodynamics, whereas the average time during which a meteoroid ablates is about 100 ms, so that the emergence of relatively volatile constituents may be diffusion-limited.
In this paper we describe a new Chemical Ablation Model (CAMOD), which contains the following processes: sputtering by inelastic collisions with air molecules before the meteoroid melts; evaporation of atoms and oxides from the molten particle; diffusion-controlled migration of the volatile constituents (Na and K) through the molten particle; and impact ionization of the ablated fragments by hyperthermal collisions with air molecules.The model is then used to calculate the injection rates of the more abundant elements as a function of altitude, to examine whether this accounts for the atomic Ca layer being depleted by more than two orders of magnitude compared with the Na layer (Plane, 2003).Finally, we will explore the curious finding that multiple metal resonance lidar observations of the same section of an individual meteor trail rarely observe more than one metal (von Zahn et al., 1999(von Zahn et al., , 2002)).

Model description
The architecture of CAMOD is illustrated with the flowchart in Fig. 1.

Thermal ablation
The classical theory of the interaction of meteoroids with the atmosphere has been described by a number of authors (e.g.Opik, 1958;McKinley, 1961;Evans, 1966;Hunten et al., 1980;Bronshten, 1983;Kalashnikova et al., 2000).The meteoroid is treated as a homogenous spherical particle.Radial heat transfer is assumed fast enough so the particle is isothermal along its whole path through the atmosphere (we explore this assumption below).The interaction with the atmosphere occurs in the free molecular flow regime in which the molecular collision mean free path is larger than the dimension of the meteoroid and thus no shock structure can develop.Then the loss of momentum of the particle of radius R and density ρ m , due to interaction with the atmosphere, is equal to the momentum of the impinging air molecules: where V is the velocity of the particle, ρ a is the mass density of the atmosphere and g the gravitational acceleration.The dimensionless "free-molecular drag" coefficient expresses the efficiency of the momentum transfer between the meteoroid and impinging air molecules.typically lies between 0.5 and 1 (Hughes, 1978).The gravitational term ρ m g is negligible in the range of meteoroid velocities.
The thermal energy received by the meteoroid from the impinging air molecules is balanced by radiative loss, temperature increase, melting, phase transitions, and by vaporization of the meteoric constituents (Jones and Kaiser, 1966).The energy conservation principle then gives: (2) The left-hand side of equation 2 represents the frictional heating term.The dimensionless parameter is the "free molecular heat transfer coefficient", which is the fraction of the incident kinetic energy which is transferred to the particle.The terms on the right-hand side of the equation represent the loss of thermal energy by the particle.The first term is the radiation loss where ε is the emissivity coefficient, σ is the Stefan-Boltzmann constant, and T and T env are the temperature of the particle surface and the atmospheric environment, respectively.The second term is the heat consumed to increase the temperature of the particle (C is the bulk specific heat).The last term is the heat consumed in the transfer of particle mass into the gas phase, where L is the latent heat of vaporization (or sublimation if the particle has not melted).The emission efficiency decreases when the size of the particle is comparable to or smaller than the radiation wavelength (Bohren and Huffman, 1983).A black body at 2000 K has a peak emission wavelength of 1.45 µm.Particles below this size contribute less than 10% to the total meteoroid mass entering the atmosphere, according to the LDEF distribution (Love and Brownlee, 1993;McBride et al., 1999), so the size dependence of the radiative heat loss term should be negligible.
T. Vondrak et al.: Meteoric ablation The change of particle height z with time is a function of its impact angle χ , defined as the angle to the zenith: The mass loss rate is calculated using Langmuir evaporation (Markova et al., 1986;Love and Brownlee, 1991;Mc-Neil et al., 1998), which assumes that the rate of evaporation into a vacuum is equal to the rate of evaporation needed to balance the rate of uptake of a species i in a closed system.The rate of mass release of species i with molecular weight µ i , from a particle of area S, is given by where the superscript A refers to thermal ablation.γ is the uptake (or sticking) coefficient, equal to the probability that the molecule is retained on the surface, or within the particle, after collision.p i is the thermodynamic equilibrium pressure of species i in the gas phase.The total mass loss rate due to ablation is then the sum over all gas-phase components: For reasons discussed below, the uptake coefficient is set to zero in the model until the temperature reaches the melting point of the particle.

Non-thermal mass loss
High velocity collisions with air molecules cause the nonthermal sputtering of atoms from the surface of the meteoroid, which results in mass loss before the particle melts and thermal evaporation dominates.We have therefore included a recent treatment of non-thermal interactions (Hill et al., 2005;Rogers et al., 2005), which was originally developed for gas-grain interactions in the interstellar medium (Tielens et al., 1994).The sputtering yields calculated in this treatment apply to interstellar ions such as H + , He + and D + , rather than the neutral constituents of interest in the MLT such as N 2 , O 2 and O.Although there are few experimental studies of sputtering from neutral bombardment, one study of low-energy neutral and ionised neon impacting on Cu, Si and SiO 2 surfaces found the sputtering yields to be similar for both the neutral and ion beams, with the neutral yield/ion yield ratio ranging from 1 for Cu to 0.7 for SiO 2 (Mizutani and Nishimatsu, 1988).Therefore, it seems reasonable to use this method in the CAMOD model.
The total sputtering yield at normal incidence Y (E, θ = 0) per projectile of energy E(in erg), due to a cascade of binary collisions in which a fraction of recoiling surface atoms acquire energy above the surface binding forces (U 0 , in eV), is given by: α is a function of the mass ratio between the target and projectile, and S n (E) is the nuclear stopping cross-section (in erg cm 2 ) which can be expressed as a function of the atomic mass and atomic number of the projectile (M p , Z p ), and the average atomic mass and atomic number of the target (M t , Z t ) (Sigmund, 1981): where e is the electron charge (1 esu), and the universal function s n (ε pt ) can be approximated as (Matsunami et al., 1980): The parameter α depends on the mass ratio of target and projectile, and in the range 0.5<M t /M p <10 can be taken as: For a mass ratio <0.5, α remains approximately constant with a value around 0.2.For a mass ratio exceeding ∼5, Eq. ( 8) overestimates the multiple scattering processes.To correct for this, α is reduced by multiplying by the ratio of the mean projected range R proj to the mean penetrated path length: K depends on the material, and a value of 0.1 has been used for an iron silicate surface (Tielens et al., 1994).The screening length of the Coulomb interaction between the nuclei, a, is a function of their atomic numbers: where a 0 is the Bohr radius (in cm).The dimensionless quantity ε pt is a function of the atomic masses and numbers, screening length and the projectile kinetic energy: Sputtering can only occur if the kinetic energy of the impinging molecule overcomes the surface binding energy, U 0 .Near the threshold the ejection of the target surface atoms is due to the recoiling projectile atoms.The threshold projectile energy is then (Bohdansky et al., 1980;Anderson and Bay, 1981): and The parameter G is equal to: However, over the range of meteoroid masses considered in the model (5×10 −18 -5×10 −3 g), only Eq. (13a) applies.When the solution of Eq. ( 6) is restricted to low projectile energies (E 20E th ), where a large collision cascade is absent and only a smaller fraction of recoiling atoms has kinetic energy above U 0 , the momentum transfer is no longer isotropic.This effect is accounted for by substituting S n for S n in Eq. ( 6) (Bohdansky et al., 1980;Yamamura et al., 1983;Bohdansky, 1984): Based on experimental evidence the angle-averaged sputtering yield for low energy impact when the cos(θ) −1 angular distribution rule holds can be taken as ∼2Y (E,θ = 0) (Draine and Salpeter, 1979).Thus for a spherical target the mass loss rate averaged over the interval 0 to π/2 is: where the superscript s refers to sputtering, the index i runs over all atmospheric components, ρ i is the number density of the i-th atmospheric component, and Y i (E,θ = 0) is the sputtering cross-section for the i-th atmospheric component.
The atmospheric components treated in the model are O 2 , N 2 , He, Ar, N, O and H.The sputtering cross section was assumed equal for all elements of which the target meteoroid was composed, and the elemental loss rate was derived using the bulk element atomic fractions (x i ) in the meteoroid: The total mass loss rate is then the sum of the thermal ablation and sputtering rates:

Ionization of the ablated fragments
The ionization coefficient of an atom (or molecule) as a function of velocity, β(V ), is calculated from the expression (Jones, 1997)  (Cuderman, 1972).The line through the points is a least-squares fit 0.347 (Bydin and Bukhteev, 1960).
The dashed line is a least squares fit 0.933 where The threshold velocity V 0 is given by where M e and φ are the mass and ionization potential of the atom (or molecule), respectively, and M a is the average molecular mass of air.Jones does not give a value of the parameter c for Ca (Jones, 1997).We have therefore calculated β 0 for Ca using the value of c for Mg (9.29×10 −6 ; Jones, 1997), but substituting the Ca ionization potential (6.113 eV; Lide, 1993).The ionization coefficients of K and Na were derived by extrapolating measurements of the collisional ionization cross sections of fast alkali metal beams (Bydin and Bukhteev, 1960;Cuderman, 1972).The β(V ) dependence was constructed as a least-squares fit of the experimental cross-sections to the expression: The relative abundances of sodium and potassium in chondritic meteorites are 0.80 and 0.051 atomic %, respectively (Mason, 1971;Fegley and Cameron, 1987;Lodders and Fegley, 1988;Sears and Dodd, 1988).Therefore, the evaporation of these elements from the meteoroid does not change the size of the particle nor affect its temperature significantly.Thus the temperature of the particle derived from Eq. ( 2) can be used to estimate the diffusion coefficient, D(T ), of the alkali metals in the molten meteoroid.The temperature dependence of D(T ) is given approximately (ignoring the activation volume) by: where D 0 is 1.48×10 −5 m 2 s −1 for K and 9.84×10 −5 m 2 s −1 for Na (Sigurdsson et al., 2000).The activation energy (E A ) for K does not seem to have been reported, and so we have adopted for K the value for Na (167 kJ mol −1 ).The fraction of a species diffusing from a spherical particle in the time interval t is given by (Crank, 1975): Equation ( 24) was solved with a vertical resolution of 100 m; t is then the time the particle travelled each 100 m step, within which T was held constant.

Vapour -molten meteoroid thermodynamic equilibrium
The vapour pressures of the species evaporating from the meteoroid were calculated using the MAGMA chemical equilibrium code (Fegley and Cameron, 1987;Schaefer and Fegley, 2004;Schaefer and Fegley, 2005).The code is robust and fast, and exhibits good agreement with experimental data for a wide range of temperatures and silicate melt compositions.MAGMA uses the ideal mixing of complex components modelled by Hastie and coworkers (Hastie et al., 1982a(Hastie et al., , b, 1984;;Hastie andBonnell, 1985, 1986), and takes account of the equilibria of eight metal oxides: SiO 2 , MgO, FeO, Al 2 O 3 , TiO 2 , CaO, Na 2 O and K 2 O.The chemical equilibria in the melts are modelled using thermodynamic activities, a i , which are given by where x i is the mole fraction of oxide in the melt and γ i is the Raoultian activity coefficient of the oxide relative to the pure liquid oxide.The melt is then modelled as a nonideal solution in which the concentrations of the unbound metal oxides are reduced by the formation of complex oxides and pseudo-components.The activities of the oxides in the non-ideal melts are calculated from the mole fractions of the unbound oxides in the melt, i.e.
where x * oxide is the mole fraction of the unbound oxide in the melt.The equilibria between components in the melt, and between the melt and vapour, are calculated simultaneously.It is assumed that the oxides evaporate stoichiometrically.The vapour pressure of each gas-phase species is then used to determine the mass loss rate of that component from the melt by invoking Langmuir evaporation (Eq.4).This approach is similar to the Equilibrium Reference Model for evaporation into a vacuum from a silicate melt that was described by Alexander (2001).It is important to note that this model reproduced the measured elemental evaporation rates from chondritic-type melts over a temperature range of 1900-2300 K, for up to 95% of mass loss.

Integration and the initial distribution of the meteoroids
A fourth-order Runge-Kutta numerical integration scheme with adaptive step size control (Press et al., 1992) is used to integrate the above sets of differential equations (Eqs.1-5, and 16-18).As shown in Fig. 1, the integration starts at an altitude of 500 km and meteoroid temperature of 200 K.After the particle has started to heat, the integration is stopped if the meteoroid mass decreases below 10 −26 g or the meteoroid temperature falls below 400 K. Otherwise the integration continues until the final altitude of 40 km.The resulting masses, ablation rates and temperature, which are calculated with a variable height resolution (maximum 22 m), are then interpolated onto a regularly spaced 100 m step altitude grid.
The model input parameters are listed in Table 1.The meteoroid mass and velocity distributions from the LDEF satellite experiment (Love and Brownlee, 1993;McBride et al., 1999) were used.The mass distribution in the model is from 5×10 −18 g to 5×10 −3 g.Larger particles would not fulfil the isothermicity condition (see Sect. 3.1 below).However, since particles heavier than 5×10 −3 g contribute less than 9% to the total input mass, imposing this upper limit will not significantly affect the calculated total mass released into the Earth's atmosphere.The extraterrestrial material is assumed to have a CI chondrite composition, with elemental abundances in Table 2 that are taken from three sources (Mason, 1971;Lodders and Fegley, 1988;Sears and Dodd, 1988).

Results and discussion
All the results presented here are for the conditions of March at 40 • N. The number densities of atmospheric constituents were calculated using the MSISE-90 model for 2001 (Hedin et al., 1991).In the following discussion, the "most probable" meteoroid refers to a particle of mass 5 µg (the most probable size in the LDEF size distribution, McBride et al., Free molecular drag coefficient, 1.0 (Hughes, 1978) Meteoroid density, ρ m 0.3, 1.0, 2.0, 3.0, or 3.5 x10 3 kg m −3 Molecular heat transfer coefficient, 1.0 Emissivity,ε 1.0 Temperature of environment, T env 200 K Specific heat of meteoroid,C 1×10 3 J kg −1 K −1 Latent heat of vaporisation,L (3.0 or 6.0)×10 6 J kg −1 Impact angle, χ 37 • or 0 • Uptake coefficient, γ 1.0 * Surface binding energy, U 0 5.7 eV (Tielens et al., 1994) Material constant (silicate), K 0.1 (Tielens et al., 1994) *The value of sticking coefficient is set to 0 for temperatures below 1700 K, 0.5 for 1800 K and 1 for temperatures higher than 1900 K. Within this interval the values follow a sigmoid curve (see Fig. 4a).

Heat transport in the meteoric particle
The application of the thermodynamic equilibrium model requires a uniform temperature across the meteoroid.If the particle is not isothermal the composition and chemical states of elements would differ between the skin and inner region of the particle.In these circumstances the MAGMA code would not provide the correct vapour pressures of the gasphase species from the bulk composition of the particle.The isothermal condition can be assessed by the dimensionless Biot number.For a spherical object the Biot number is given by where h is the heat transfer coefficient (W m −2 K −1 ), R is the radius and k the thermal conductivity (W m −1 K −1 ).If Bi<0.1, the object can be considered isothermal (Kakac and Yener, 1985).If non-radiative heat losses from the particle are omitted, the heat transfer coefficient can be expressed as (Love and Brownlee, 1991): Figure 3 shows the calculated Bi for a range of initial particle masses over the temperature range in which the meteoroid ablates, using a thermal conductivity for olivine rock of 1 W m −1 K 1 (Buttner et al., 1998).This shows that the most probable 5 µg particle is isothermal up to a very high temperature, and the largest mass used in the model stays isothermal up to ∼2000 K if the density is 2 g cm −3 .Because meteoroids lose mass extensively above ∼2800 K they become much smaller at the highest temperatures attained, reducing Bi further.Thus the isothermicity condition seems fulfilled for the whole mass distribution up to 5×10 −3 g.Table 2. Initial composition of meteoroids (Mason, 1971;Fegley and Cameron, 1987;Lodders and Fegley, 1988;Sears and Dodd, 1988) a Elemental abundance on the cosmochemical scale: Si is taken as 1.00×10 6 .

Melting of the meteoroid
The MAGMA code gives non-zero vapour pressures for metal species at temperatures below the melting temperature of the meteoroid.However, thermodynamic equilibrium in the meteoroid cannot be established before the particle melts.
Thus the melting temperature has to be entered as an external parameter into the integration procedure.Mineralogically, CI chondrites have an olivine composition (Mason, 1971;Lodders and Fegley, 1988;Sears and Dodd, 1988).The olivine phase diagram shows that for the chondritic Fe/Mg ratio of ∼0.8, olivine starts to melt at 1730 K (Fig. 4a).The onset of particle evaporation in CAMOD is simulated by applying a sigmoid temperature dependence to the uptake coefficient γ in Eq. ( 4).A finite smooth melting range describes in a more realistic way the phase change than an instantaneous transition to the molten state at the threshold temperature, and also prevents instability in the numerical integration.In The dependence of the isothermicity condition on the ablation temperatures, for a range of particle masses (0.5 to 500 µg) and particle densities (solid lines 1 g cm −3 , broken lines 3 g cm −3 ).For the points below the curves the Biot number is lower than 0.1 and thus the isothermicity condition is fulfilled.The dotted line marks the thermal conductivity of olivine.
practice, the value of the melting temperature only affects the altitude of the deposition of the alkali metals: within the temperature range which encompasses the olivine phase equilibrium (1600-2000 K), the altitude of the Na ablation peak varies from 96-101 km (Fig. 4b).We have adopted a meteoroid melting temperature of 1800 K, and γ varies from 0 (T <1700 K) to 1 (T >1900 K) for all the model runs discussed below.The meteoroid mass-velocity combinations which undergo thermal ablation are shown in Fig. 5, for a range of particle densities (0.3-3 g cm −3 ).This demonstrates that the velocity at which the particle temperature reaches the ablation threshold is a steep function of particle density.
Figure 6 shows the dependence with entry velocity of the height at which meteoroids (mass range 0.05-500 µg) are predicted to melt.The modelled dependence is compared with the appearance height of meteor head echoes observed by the 430 MHz Arecibo incoherent radar (Janches and ReVelle, 2005).Note that the velocity measured by the radar is the vertical component, along the line-of-sight of the radar.There is very good agreement for velocities between 20 km s −1 and 45 km s −1 .At velocities below 20 km s −1 the head echoes appear 10-15 km higher.However, this discrepancy is most probably explained by meteoroids crossing the radar beam at large entry angles, and thus with small measured vertical velocity components.Indeed, this must explain the velocities below 11 km s −1 observed by the radar.At velocities above 45 km s −1 , the model predicts much higher appearance heights for smaller meteoroids than the average appearance heights of the Arecibo observations.There are two likely reasons for this.First, immediately after the particle melts, the volatile elements Na and K ablate (see Figs. 8 and 10 below).Although these elements ionize efficiently at high impact velocities (see Fig. 14 below), they comprise less than 1% of the meteoroid mass and so the rate of electron production contributing to the meteor head echo is small.In contrast, the bulk constituents (Fe, Mg and Si) ablate about 10 km lower (Fig. 10, below), when the rate of electron production will be considerably higher.Second, the head echo is proportional to the square root of the electron density around the meteoroid, and the electron density in turn varies inversely with atmospheric density because the electron cloud expands to about one mean free path around the meteoroid (Dyrud et al., 2005).The head echo thus drops by Fig. 5. Dependence of the thermal ablation threshold on meteoroid mass and velocity for a range of particle densities (0.3-3 g cm −3 ).Only particles to the right of the curves reach a temperature greater than 1800 K.
∼10 dB for a 10 km height increase.Note that even though these two factors mitigate against observing small high velocity meteoroids at the altitudes where thermal ablation begins, the Arecibo radar does in fact see a small number of particles above 120 km (Janches et al., 2003).The coupling of CAMOD with a meteor head echo model will be used for exploring these effects in a future study.

Non-thermal loss
Sputtering as parallel mechanism to ablation was first proposed by Öpik (1958), and the main mass loss mechanism for particles which do not melt (Lebedinets and Shushkova, 1970a, b).Sputtering has been proposed as an explanation for the luminosity of high-altitude meteors above 130 km , which cannot be explained by classical ablation theory (Hill et al., 2004;Popova et al., 2007;Vinkovic, 2007).Several other processes have been observed to occur during collisions between high-energy beams and various surfaces in laboratory experiments (Behrisch and Eckstein, 2007), such as neutralisation of ionic projectiles followed by their reflection (Kleyn, 1992), negative ion formation (Greber, 1997), and compaction of the impacted surface (Mizutani, 1995).However, only sputtering is considered in the present study since we are primarily concerned with mass loss from a meteoroid.In 2002, a simple sputtering mechanism was incorporated into a meteor model (Coulson, 2002;Coulson and Wickramasinghe, 2003), followed by a more detailed sputtering model which included the variation μ μ μ μ μ Fig. 6.Comparison of the average appearance heights of meteor head echoes observed by the 430 MHz radar at Arecibo as a function of entry velocity (points, with 1σ standard deviations), and the modelled appearance heights (lines) defined as the altitude when the particle temperature reaches 1800 K and the alkali metals start to ablate thermally, for a range of particle masses (0.05-500 µg).The recorded velocities lower than 10 km s −1 (hollow points) correspond to particles arriving at large entry angles, which therefore have a small vertical velocity component.
of atmosphere composition along the meteoroid trajectory (Rogers et al., 2005).This later model used the treatment developed for ion impact on surfaces by Tielens et al. (1994), which we have also incorporated into CAMOD.Notwithstanding the fact that this treatment was developed for ionic projectiles rather than the neutral projectiles that predominate in the thermosphere (Section 2.2), the results discussed below should provide a useful indication of the relative importance of sputtering.
The fraction of mass loss before the meteoroid reaches 1700 K and thermal ablation starts is shown in Fig. 7, over the range of meteoroid masses and entry velocities.Clearly the effect of sputtering is negligible for meteoroids with velocities up to 45 km s −1 and heavier than 10 −11 g.Only the fastest particles up to masses of 10 −14 g are reduced to approximately half of their initial mass by sputtering.

Diffusion-controlled ablation of Na and K
The alkali elements are the most volatile metallic constituents of meteoroids and, according to the thermodynamic model with Langmuir evaporation, should ablate very quickly after the meteoroid melts.This raises the question of whether the rate of ablation of Na (and K)   be kinetically controlled by diffusion out of the molten particle.Figure 8 shows the altitude profiles of Na ablating from the most probable meteoroid for a range of entry velocities, when evaporation is controlled either by diffusion (Eqs.23 and 24) or by Langmuir evaporation (Eq.4).Inspection of this figure demonstrates that for slow meteoroids (V up to ∼20 km s 1 ) the ablation profiles are similar: the FWHMs (full-width-at-half-maximum) of the ablation profiles differ by less than 50%.In contrast, radial diffusion significantly broadens (by a factor of ∼2.5) the ablation profile of the fastest meteoroid at 72 km s −1 .However, the weighting of the fast particles in the meteoroid velocity distribution (Love and Brownlee, 1993;McBride et al., 1999) is a few orders of magnitude lower than the weighting of the most probable V (∼20 km s 1 ).Thus, the total injection profiles of the alkali metals, obtained by summing over the whole meteoroid velocity and mass distribution (Sect.3.5), is hardly affected by including diffusion control.

Elemental ablation profiles in the atmosphere
The fraction of a meteoroid that ablates is shown as a function of mass and entry velocity in Fig. 9.The lower mass threshold for complete evaporation of fastest meteoroid particles is ∼10 −10 g.The mass threshold for nearly complete (95%) ablation as a function of meteoroid velocity can be expressed by the following inequality (ρ m = 2 g cm −3 ): where V is the meteoroid velocity in km s −1 .The elemental deposition profiles for the most probable meteoroid are shown in Fig. 10.Na and K are released at ∼99 km; their vaporization is almost complete before other elements start to evaporate.Indeed, the less volatile major components Si, Fe and Mg are released 10-15 km lower in the atmosphere.The ablation of Na, K, Fe, Si, and Mg is nearly complete and thus the relative abundance of these elements in the gas phase is very similar to the meteoroid before entry.The behaviour of Ca is strikingly different: less than 1% of this element is released from the particle (Table 3).The Na/Ca ratio in the gas phase is thus ∼170 times larger than in the meteoroid.Interestingly, this ratio is similar to the ratio of the column densities of the background Na and Ca layer observed by lidar (Qian and Gardner, 1995).
In contrast, the integration of the elemental ablation profiles over the meteoroid velocity and mass distributions, to yield the total injection rate of each element, presents a somewhat different picture.Figure 11 shows that the peak release of the alkali metals remains approximately 10 km above the major meteoric components, but the Na/Ca ratio is now only 4:1 (Table 4).This arises because the fraction of Ca that ablates is strongly dependent on entry velocity.As shown in Fig. 12 (middle panel), Ca is not released from the most probable meteoroid until V exceeds 20 km s −1 .At higher velocities, the fraction of ablated Ca increases steeply to reach 1.0 at V = 35 km s −1 .Fig. 11.Injection rates of individual elements, integrated over the LDEF distributions of meteoroid mass (5×10 −18 -5×10 −3 g) and velocity (11.5-72.5 km s 1 ).Note that these injection rates should be scaled by a factor of 0.34 if the total mass flux is 30 t d −1 (see text).
The threshold velocity for Ca ablation and the velocity at which complete ablation of Ca occurs also depend strongly on the meteoroid mass, as shown in the three panels of Fig. 12.A 50 µg particle has an ablation threshold velocity of only 15 km s −1 and the Ca ablates completely above 30 km s −1 , whereas for a 0.05 µg particle these velocities are 35 and 55 km s −1 , respectively.Indeed, simultaneous common-volume lidar measurements of Ca and Fe have shown that the Ca/Fe ratio in fresh meteor trails approaches the chondritic ratio as the meteoroids get larger (von Zahn et al., 2002).These model predictions are also in accord with steady-state evaporation experiments using FeO-MgO-SiO 2 -CaO-Al 2 O 3 chondritic melts (Hashimoto et al., 1979;Hashimoto, 1983).The first stage of heating is associated with the loss of FeO, and a significant enrichment of CaO and Al 2 O 3 occurs after 70% of the mass has evaporated.Lastly, cosmic spherules that exhibit progressive heating are characterized by loss of the more volatile oxides and become Ca-Al-Ti-rich (Taylor et al., 2000).
The injection rate profiles plotted in Fig. 11 are the current best estimates for use in atmospheric models of the mesospheric metal layers.Note that these injection rates should be scaled by a factor of 0.34 if the total mass flux is around 30 t d −1 , as argued in Sect. 1. Future refinements to this average picture can be made by applying CAMOD to the seasonal and latitudinal variations of the meteoroid mass, velocity and entry angle distributions, as well as considering different compositions for meteoroids depending on their origin (Janches et al., 2006;Fentzke and Janches, 2008).3.6 Lidar detection of K, Ca and Ca + in meteor trails One test of CAMOD is to explain the finding that when common-volume lidar measurements are made of two or three metals, meteor trails that pass through the common volume rarely contain more than a single metal: two metals are seen <1% of the time, and three metals only once out of thousands of observed trails (von Zahn et al., 2002).
Here we focus on simultaneous common-volume measurements of K and Ca performed during a meteor shower et al., 1999).Because lidar were zenith-pointing, and meteor trails are typically inclined at an angle the zenith, lidars would only have at most a 100 m along the trail.The most striking feature of the observations of ∼30 meteor trails was an absence of a simultaneous detection of the two elements, and a pre-Table 3. Ablated fraction and the ablation profile parameters (centroid height and FWHM) for ablation of the major elements from the most probable meteoroid (mass=5µg, entry velocity=20.5 km s −1 , χ =37 • , ρ m =2.0 g cm −3 ).ponderance (>80%) of Ca-containing trails within a height range of 92-95 km .CAMOD predicts that a 50 µg meteoroid (i.e.large enough to release a measurable amount of Ca) will inject Ca at these altitudes if the entry velocity is ∼50 km s −1 (Fig. 13, left-hand panel).In contrast, K will then ablate at an altitude of ∼112 km if Langmuir evaporation is rate-limiting, or ∼109 km when diffusion out of the meteoroid is rate-limiting.However, above 100 km the probability of detecting K in the trail is strongly reduced because of the rapid expansion of the trail by radial diffusion at the lower atmospheric pressure (Höffner et al., 1999).Furthermore, Figure 14 shows that most of the potassium (∼80%) will be ionized by hyperthermal collisions at this entry velocity, whereas a much smaller fraction (∼25%) of calcium will be ionized.To release K around 90 km, where K trails have been observed (Gerding et al., 1999), the initial velocity of the meteoroid must be around 12 km s −1 (Fig. 13, right-hand panel).At this low velocity only ∼10% of the ablating K will undergo collisional ionization (Fig. 14).However, the peak temperature reached by the meteoroid is only ∼2300 K, and so only a trace (<0.1%) of Ca ablates (Fig. 13, right panel).Thus, CAMOD is able to explain the presence of either K or Ca, but not both metals, in simultaneous lidar observations of the same section of a meteor trail.
von Zahn et al. (2002) have reported simultaneous common-volume measurements of Ca and Ca + in meteor trails, using the resonance lidar technique (Alpers et al., 1996).The mean height at which these layers were observed was 94.5 km , and the Ca + /Ca ratio was greater than 1.As shown in Fig. 14, a 50% ionization rate of Ca is only reached for impact velocities of ∼70 km s −1 , and it is very unlikely that all the trails reported by von Zahn et al. (2002) were produced by such high speed meteoroids (indeed, the lefthand panel of Fig. 13 indicates that an entry velocity around 50 km s −1 should produce a Ca trail at 95 km).This discrepancy probably indicates that our estimate of the collisional ionization efficiency of Ca, made using a semi-empirical formalism (Eqs. 19-21;Jones, 1997), underestimates β by a factor of 2 -3.
www.atmos-chem-phys.net/8/7015/2008/Fig. 15.Dependence of elemental ablation on the particle density.The ablated fraction of the initial mass (upper panel) and centroid of the ablation profile (lower panel) are shown for selected elements, integrated over the LDEF distributions of mass (5×10 −18 -5×10 −3 g) and velocity (11.5-72.5 km s −1 ).The dashed vertical line indicates the particle density inferred from the analysis of impact craters in the LDEF experiment.

Differential ablation of Ca: sensitivity to model parameters
The model results presented in Tables 3 and 4 were obtained assuming a bulk meteoroid density of 2 g cm −3 and heat of vaporization of 6 MJ kg −1 .Although a density of 2 g cm −3 was assumed for analysis of the LDEF cratering experiment (Love et al., 1995), the bulk porosity of incoming extraterrestrial particles is rather uncertain (Love et al., 1993;Love et al., 1994;Rietmeijer, 2002).We now consider a series of densities in the range 0.3 to 3.5 g cm −3 .The lowest density is typical of highly porous extraterrestrial material of cometary origin, whereas asteroidal material is thought to have a density of ∼3 g cm −3 (Rietmeijer, 2002).For comparison, the density of olivine lies between 3.2 g cm −3 (forsterite) and 4.3 g cm −3 (fayalite) (Lide, 1993).
Figure 15 shows the effect of the meteoroid density on the ablated fraction of each element and the centroid height of its ablation profile, integrated over the LDEF particle mass and velocity distributions.When the density increases from 0.3 to 3.5 g cm −3 (the mass distribution is conserved, so the particle volume decreases as ρ −1 m ), the total ablated fraction increases nearly 2-fold.However, the increases of the individual elements vary from 1.3 (Na, K) to 2.4 (Ca).The centroid heights of the elemental ablation profiles decrease by 3.6 km (Na, K) to 7.1 km (Ca).Thus the relative amount of ablated Ca changes by less than a factor of 2 for a 10-fold increase of ρ m .This indicates that a meteoroid population with a density very different from 2 g cm −3 is unlikely to explain the depletion of atomic Ca, relative to Na and Fe, in the upper mesosphere.
Vaporization of the meteoroid is the main cooling mechanism once the particle has melted.In spite of this, the effect of halving the heat of vaporization, L, brings about a less than 5% increase of the deposited metals and a negligible shift of the injection profile peaks (Table 5).The ablation efficiency is also not sensitive to the specific heat of the meteoroid, C. For example, increasing C by a factor of 3 causes less than 1% change in the total ablated fraction, and shifts the ablation centroid heights less than 1 km .
Lastly we consider the meteoroid incidence angleχ.Table 4 shows that the injection rates of the major components and alkali metals are only slightly larger (<10%) for zenith incidence compared to χ=37 • , although the refractive metals (Ca, Ti, and Al) show more sensitivity toχ: their integrated injection rates increase by 13-20%.The centroid heights of the injection profiles shift down between 0.1 and 3 km at χ = 0 • .

Conclusions
A model of the chemical ablation of meteoroids in the upper atmosphere has been described.Na and K evaporate first, several km above the main constituents Fe, Mg and Si.CAMOD successfully predicts over a broad of meteoroid velocities the meteor head echo appearance heights measured by incoherent scatter radar, although the discrepancy noted at the highest velocities will be the subject of future work where CAMOD is coupled to a head echo model.CAMOD also explains common-volume lidar observations of K, Ca and Ca + in fresh meteor trails.The ablation of Ca relative to the major constituents and alkalis depends strongly on the initial velocity of meteoroids.The Na:Ca ratio of 100-360:1 observed in the metal atom layers in the upper mesosphere is reproduced only for slow meteoroids (velocity<20 km s −1 ).However, CAMOD predicts a ratio of only 4:1 integrated over the commonly assumed distribution of meteoroid mass and velocity.This result is also not very sensitive to the meteoroid density, heat of evaporation, and angle of entry assumed in the model.The explanation for the more than 2 orders of magnitude depletion of Ca relative to Na atoms in the mesosphere, compared with their chondritic abundance (∼1:1), must therefore lie elsewhere.

Fig. 1 .
Fig. 1.Flow chart illustrating the architecture of the chemical ablation model.
22) where A = 0.933 and B =-1.94 for Na, and A = 0.347 and B =-1.702 for K.The extrapolated ionization coefficients and experimental values for the alkali metals are shown in Fig. 2 together with the experimental values.The β values for other metals derived from Eq. (19) are included for comparison.T. Vondrak et al.: Meteoric ablation 2.4 Diffusion-controlled ablation of alkali metals

T.Fig. 3 .
Fig.3.Heat transport in the meteoroid.The dependence of the isothermicity condition on the ablation temperatures, for a range of particle masses (0.5 to 500 µg) and particle densities (solid lines 1 g cm −3 , broken lines 3 g cm −3 ).For the points below the curves the Biot number is lower than 0.1 and thus the isothermicity condition is fulfilled.The dotted line marks the thermal conductivity of olivine.

Fig. 4 .
Fig. 4. Melting temperature of the meteoroid.a: Phase diagram of olivine and the assumed melting temperature of the chondritic particle.The dashed line represents the temperature dependence of the uptake coefficient γ , which is used to implement a smooth melting transition in the model.b: The dependence of the Na ablation profile (solid lines) on the meteoroid melting temperature, for a 5 µg meteoroid at 20 km s −1 .The Fe ablation profile (broken line) for 1800 K is shown for comparison.

Fig. 7 .Fig. 8 .
Fig.7.Non-thermal mass loss.The fraction of a particle sputtered before it reaches 1700 K, as a function of meteoroid mass and velocity.Meteoroids smaller than ∼10 −13 g do not reach the melting temperature (1800 K), so the sputtered fraction is the whole mass loss.Particle density = 2 g cm −3 .

Fig. 9 .
Fig.9.Ablated fraction of a meteoroid (particle density 2 g cm −3 ), as a function of mass and velocity (the wavy structure of the contours is an artefact caused by the discretisation of the particle mass distribution).

Fig. 10 .
Fig. 10.Ablation profiles of individual elements from a 5 µg meteoroid entering at 20 km s −1 .The particle temperature is shown with the solid line, referenced to the top abscissa.

Fig. 14 .
Fig. 14.Dependence of the ionized fraction β of selected elements (at their peak ablation altitude) on the initial velocity of a 5 µg meteoroid.

Table 4 .
Ablated fraction and the ablation profile parameters (centroid height and FWHM) for ablation of the major elements from meteoroids integrated over the whole mass and velocity distribution.

Table 5 .
Dependence of the elemental ablated fractions and centroid heights on the heat of vaporisation (L) for a meteoroid density 2 g cm −3 .