Silicon chemistry in the mesosphere and lower thermosphere.

Silicon is one of the most abundant elements in cosmic dust, and meteoric ablation injects a significant amount of Si into the atmosphere above 80 km. In this study, a new model for silicon chemistry in the mesosphere/lower thermosphere is described, based on recent laboratory kinetic studies of Si, SiO, SiO2, and Si+. Electronic structure calculations and statistical rate theory are used to show that the likely fate of SiO2 is a two-step hydration to silicic acid (Si(OH)4), which then polymerizes with metal oxides and hydroxides to form meteoric smoke particles. This chemistry is then incorporated into a whole atmosphere chemistry-climate model. The vertical profiles of Si+ and the Si+/Fe+ ratio are shown to be in good agreement with rocket-borne mass spectrometric measurements between 90 and 110 km. Si+ has consistently been observed to be the major meteoric ion around 110 km; this implies that the relative injection rate of Si from meteoric ablation, compared to metals such as Fe and Mg, is significantly larger than expected based on their relative chondritic abundances. Finally, the global abundances of SiO and Si(OH)4 show clear evidence of the seasonal meteoric input function, which is much less pronounced in the case of other meteoric species.

Introduction

Silicon is one of the most abundant elements in cosmic dust, comprising around 11% by mass (9% elemental abundance) [Asplund et al., 2009]. The amount of cosmic dust entering the Earth's atmosphere each day is still quite uncertain [Plane, 2012], with recent estimates in the range of 14 to 150 t d−1 [Carrillo‐Sánchez et al., 2015; Huang et al., 2015]. The fraction of silicon which ablates from this dust in the upper atmosphere is also highly dependent on the dust size and entry velocity distributions [Carrillo‐Sánchez et al., 2015]. Although neutral silicon compounds of extraterrestrial origin have not yet been detected in the mesosphere/lower thermosphere (MLT) region between 70 and 120 km, Si+ ions have been measured by rocket‐borne mass spectrometry over several decades (see reviews by Kopp et al. [1995] and Grebowsky and Aikin [2002]).

Figure 1a shows vertical profiles of Si+, Fe+, and Mg+ averaged from five rocket flights conducted by Kopp and coworkers. Two of these rockets were launched from Wallops Islands (37.8°N, USA) at 11:58 LT on 12 August 1976 (payload 18.1006 [Herrmann et al., 1978]) and at 14:03 LT on 1 January 1977 (payload 18.1008 [Meister et al., 1978]). The rest were launched from Kiruna (67.8°N, Sweden) at 01:38 LT on 13 August 1978 (payload S26/2 [Kopp et al., 1985]), at 01:44 LT on 30 November 1980 (payload 33.009 [Kopp, 1997]), and at 01:32 LT on 3 August 1982 (payload S37/P [Kopp et al., 1984]). Care needs to be taken when assigning mass‐to‐charge ratio (m/z) 28 to Si+, as N2 + produced by energetic particle precipitation in the lower thermosphere, or H2CN+ from the mesosphere, have the same mass‐to‐charge ratio. Here we have confirmed that the 29Si/28Si isotopic ratio for each of the profiles that were selected for Figure 1 was in agreement with the natural abundance ratio of 0.051 over the altitude range 90–110 km [Zbinden et al., 1975].

Figure 1

Rocket‐borne mass spectrometric observations of Si+, Fe+, and Mg+ (see text for flight details): (a) averaged vertical profiles of five rocket flights; (b and c) averaged Si+/Fe+ and Mg+/Fe+ ratios versus altitude (squares), along with the maximum and minimum ratio values (dashed lines). The thick solid lines show the corresponding WACCM modeled ratios, with minimum and maximum values indicated by thin solid lines. The dash‐dot‐dot vertical lines indicate chondritic ratios.

Inspection of Figure 1 shows that above 100 km the concentrations of Si+, Fe+, and Mg+ are similar, as expected because they have similar cosmic abundances [Asplund et al., 2009] and should ablate with similar efficiencies [Carrillo‐Sánchez et al., 2015]. However, whereas the Fe+ and Mg+ profiles maintain a fairly constant ratio over the entire height range, the relative Si+ concentration decreases by about 2 orders of magnitude below 97 km. This implies that Si+ is neutralized more easily than other meteoric metal ions in the upper mesosphere, which is an important test of an atmospheric model. The resulting neutral silicon oxides probably then polymerize with Fe and Mg oxides/hydroxides to form meteoric smoke particles (MSPs), a process which has been shown to occur rapidly in the laboratory [Saunders and Plane, 2011]. There is tentative spectroscopic evidence, obtained from the Solar Occultation For Ice Experiment (SOFIE) instrument on board the Aeronomy of Ice in the Mesosphere (AIM) satellite, that MSPs below 85 km are composed of Fe‐Mg‐silicates [Hervig et al., 2012]. Lastly, rocket‐borne measurements of negative ions below 90 km have identified an ion of mass 76 which could be either CO4 − or SiO3 −, though there are strong arguments for the latter [Viggiano et al., 1982].

The purpose of this paper is to develop a comprehensive model of the ion and neutral chemistry of silicon. Following a series of laboratory kinetic studies on the reactions of Si, SiO, and Si+ [Eyet et al., 2010; Gómez Martin et al., 2009a; Gómez Martin et al., 2009b; Gómez Martín and Plane, 2011], it is now feasible to make a significant advance on the original silicon model of Fritzenwallner and Kopp [1998]. Here we first demonstrate using high‐level theory that silicic acid (Si(OH)4) is a likely sink for silicon in the MLT. The rate coefficients for a small number of reactions, not hitherto measured, are then calculated using theoretical techniques, before the silicon chemistry scheme is placed into the Whole Atmosphere Community Climate Model (WACCM) [Marsh et al., 2007] with a meteoric input function of Si. We have used this model previously to study iron and magnesium chemistry in the MLT [Feng et al., 2013; Langowski et al., 2015]. The modeled Si+, Fe+, and Mg+ ratios are shown to compare well against rocket measurements, before the model is used to predict the speciation of neutral silicon in the MLT.

Chemistry of Silicon in the MLT

Figure 2 is a schematic diagram of the ion‐molecule and neutral chemistry of silicon in the MLT. The rate coefficients of the reactions shown with black arrows have been measured in the laboratory (see footnotes to Table 1); rate coefficients for the dissociative recombination reactions of silicon‐containing ions with electrons (dashed grey arrows) have been set to a typical value for this class of reactions [Florescu‐Mitchell and Mitchell, 2006]; and the rate coefficients for the reactions shown with red arrows have been calculated, as described below.

Figure 2

Schematic diagram of the important reaction pathways for silicon produced by meteoric ablation in the MLT. Neutral and ionic Si species are enclosed in green and blue boxes, respectively. Solid black arrows depict chemical reactions which have been studied in the laboratory; dashed line arrows indicate dissociative electron recombination reactions (see text); solid red arrows depict reactions where the rate coefficients are calculated theoretically in this study. The main neutral sink is Si(OH)4, which is an important constituent of meteoric smoke (blue arrow).

Table 1

Neutral and Ionic Gas Phase Reactions in the Silicon Model

Reaction Number	Reaction	Rate Coefficienta	Source
		Neutral Chemistry
(R1)	Si + O2 → SiO + O	9.5 × 10−11 + 1.8 × 10−10exp(−T/115)	Gómez Martin et al. [2009a]
(R2)	Si + O3 → SiO + O2	4.0 × 10−10	Gómez Martin et al. [2009a]
(R3)	SiO + O3 → SiO2 + O2	5.0 × 10−13	Gómez Martin et al. [2009b]
(R4)	SiO + OH → SiO2 + H	log10(k) = −4.921 − 3.801 log10(T) + 0.5023(log10(T))2	Gómez Martin et al. [2009b] and Plane [2013]
(R5)	SiO2 + H → SiO + OH	log10(k) = −19.83 + 7.323 log10(T) − 1.366 (log10(T))2	Plane [2013]
(R6)	SiO2 + H2O (+ M) → OSi(OH)2	k rec,0 = 5.54 × 10−25 (T/200)−3.69	Calculated from RRKM theory (see text)
		k rec,∞ = 3 × 10−10 (T/200)1/6
		F c = 0.188
(R7)	OSi(OH)2 + H2O (+ M) → Si(OH)4	k rec,0 = 2.71 × 10−21 (T/200)−4.22	Calculated from RRKM theory (see text)
		k rec,∞ = 3 × 10−10 (T/200)1/6
		F c = 0.249
		Ionic Chemistry
(R8)	SiO + O2 + → SiO+ + O2	4.2 × 10−9 (T/200)−0.5	Calculated from Langevin theory (see text)
(R9)	SiO+ + O → Si+ + O2	6.0 × 10−10	Calculated from Langevin theory (see text)
(R10)	Si+ + O3 → SiO+ + O2	3.3 × 10−10	Gómez Martín and Plane [2011]
(R11)	Si+ + O3 → SiO + O2 +	3.2 × 10−10	Gómez Martín and Plane [2011]
(R12)	Si+ + H2O → SiOH+ + H	2.5 × 10−10	Averaged value from Fahey et al. [1981], Wlodek et al. [1987], and Glosik et al. [1995]
(R13)	Si+ + O2(1Δg) → SiO+ + O	3.6 × 10−11(T/200)−2.2	Eyet et al. [2010]
(R14)	Si+ + O2 (+M) → SiO2 +	9.0 × 10−29 (T/200)−2.34	Gómez Martín and Plane [2011], increased by a factor of 3 for N2 as the third body rather than He
(R15)	SiO+ + O3 → SiO2 + O2 +	6.0 × 10−10	Gómez Martín and Plane [2011]
(R16)	SiO2 + + O → SiO+ + O2	6.0 × 10−10	Calculated from Langevin theory (see text)
(R17)	SiO+ + H2O → SiOH+ + OH	2.9 × 10−9 (T/200)−0.5	Calculated from Langevin theory (see text)
(R18)	SiO+ + H2 → SiOH+ + H	3.2 × 10−10	Fahey et al. [1981]
(R19)	SiO2 + + H → SiOH+ + O	2.0 × 10−10	Calculated from Langevin theory (see text)
(R20)	SiO+ + e − → Si + O	3.0 × 10−7	Set to typical value for dissociative recombination reactions of small molecular ions [Florescu‐Mitchell and Mitchell, 2006]
(R21)	SiOH+ + e − → Si + OH	3.0 × 10−7	Set to typical value for dissociative recombination reactions of small molecular ions [Florescu‐Mitchell and Mitchell, 2006]
(R22)	SiO2 + + e − → Si + O2	3.0 × 10−7	Set to typical value for dissociative recombination reactions of small molecular ions [Florescu‐Mitchell and Mitchell, 2006]
(R23)	Si+ + e − → Si + hv	9.2 × 10−12 (T/200)−0.56	Badnell [2006]

Units: unimolecular, s−1; bimolecular, cm3 molecule−1 s−1; termolecular, cm6 molecule−2 s−1.

Neutral Chemistry

The Chemical Ablation Model (CABMOD) [Vondrak et al., 2008] predicts that silicon mostly ablates from meteoroids in the form of SiO and SiO2, based on the thermodynamic model MAGMA [Schaefer and Fegley, 2005]. These oxides will most likely then be dissociated to atomic Si through hyperthermal collisions with air molecules, since they have initial velocities in excess of 11 km s−1. However, Si will immediately be oxidized to SiO by reaction with O2 or O3 [Gómez Martin et al., 2009a]. In turn, SiO is oxidized more slowly to SiO2 by reaction with O3 or OH [Gómez Martin et al., 2009b]. SiO2 is a relatively stable molecule: although the reaction SiO2 + O → SiO + O2 is exothermic by 63 kJ mol−1, quantum calculations indicate the existence of a high barrier (104.7 kJ mol−1) which will make this reaction very slow at mesospheric temperatures [Gómez Martin et al., 2009b]. The reaction SiO2 + H → SiO + OH is endothermic by only 6 kJ mol−1, and detailed balance with the measured rate coefficient for SiO + OH has been used to estimate its rate coefficient [Plane, 2013], which is listed in Table 1.

Here we use quantum theory calculations to explore the hydrolysis of SiO2, in order to determine whether the sequential addition of two H2O molecules can form the stable molecule Si(OH)4 under the conditions of the MLT: where M is a third body (N2 or O2). Electronic structure calculations using the Gaussian 09 suite of programs [Frisch et al., 2009] were used to map the stationary points on the electronic potential energy surfaces (PESs) of these reactions. Molecular geometries were first optimized using hybrid density functional theory, which includes some exact Hartree‐Fock exchange. The B3LYP method was used together with the 6‐311+G(2d,p) triple zeta basis set. This large, flexible basis set has both polarization and diffuse functions added to the atoms. After the optimized geometries were checked for wave function stability, the resulting rotational constants and vibrational frequencies were used in the master equation calculations described below. More accurate energies (±15 kJ mol−1 [Foresman and Frisch, 1996]) were then determined using the Complete Basis Set (CBS‐Q) method of Petersson and coworkers [Montgomery et al., 2000]. All reaction enthalpy changes reported in this paper were calculated at this level of theory.

Figure 3a illustrates the PES for reaction (1). This shows that after initially forming a SiO2‐H2O complex, rearrangement over a barrier submerged below the energy of the entrance channel leads to the product where three O atoms are now bound to the Si. The PES for reaction (2) is similar, as shown in Figure 3b: a shallow OSi(OH)2‐H2O complex rearranges over a submerged barrier to form Si(OH)4. The rate coefficients for these reactions were then estimated using Rice‐Ramsperger‐Kassel‐Markus (RRKM) theory, employing a solution of the master equation (ME) based on the inverse Laplace transform method [Robertson et al., 2007]. The calculations were carried out with the open source program MESMER (Master Equation Solver for Multi‐well Energy Reactions) [Robertson et al., 2012], which determines the temperature‐ and pressure‐dependent rate coefficient from the full microcanonical description of the system time evolution by performing an eigenvector/eigenvalue analysis similar to that described by Bartis and Widom [1974]. The internal energies of the intermediates on the PES were divided into a contiguous set of grains (width 50 cm−1) containing a bundle of rovibrational states. Each grain was then assigned a set of microcanonical rate coefficients linking it to other intermediates, calculated by RRKM theory. For dissociation to reactants, microcanonical rate coefficients were determined using inverse Laplace transformation to link them directly to the capture rate coefficient, which was set to 3 × 10−10 (T/200 K)1/6 cm3 molecule−1 s−1, where the small positive temperature dependence is characteristic of a long‐range potential governed by dispersion forces [Georgievskii and Klippenstein, 2005].

Figure 3

(a) Potential energy surface for the recombination of SiO2 and H2O to form OSi(OH)2; (b) potential energy surface for the recombination of OSi(OH)2 and H2O to form Si(OH)4. Color scheme: Si (grey), O (red), and H (white). TS = transition state.

The probability of collisional transfer between grains was estimated using the exponential down model, where the average energy for downward transitions was set to < ΔE > down = 300 cm−1 [Gilbert and Smith, 1990]. Because reactions (1) and (2) are fast and hence well into the falloff region between third‐ and second‐order kinetics even at the low pressures of the MLT, the MESMER calculations of k rec over a range of T and p were fitted to the standard Lindemann expression containing the low‐ and high‐pressure limiting rate coefficients, k rec,0 and k rec,∞, modified by a broadening factor F c [Sander et al., 2011]. The fitted values are given in Table 1. The molecular parameters used in the calculations are listed in Table 2.

Table 2

Molecular Properties of the Stationary Points on the Potential Energy Surfaces in Figure 3, Calculated at the B3LYP/6‐311+g(2d,p) Level of Theory

Molecule	Geometry (Cartesian Coordinates in Å)	Rotational Constants (GHz)	Vibrational Frequencies (cm−1)
SiO2	Si,0,0,0	6.912	298 (x2), 991, 1438
	O,1.512,0,0
	O,−1.512,0,0
SiO2‐H2O complex	Si,−0.249,0.279,−0.005	7.890, 6.934, 3.701	176, 304, 328, 376, 394, 441, 736, 992, 1395, 1569, 3719, 3835
	O,−1.475,−0.619,0.067
	O,0.559,1.568,−0.006
	O,1.151,−1.028,−0.133
	H,2.031,−0.735,0.149
	H,0,0.933,−1.916,0.190
TS from SiO2‐H2O to OSi(OH)2	Si,−0.224,0.072,0.0169	10.747, 5.944, 3.853	1227 i, 293, 314, 486, 510, 748, 798, 966, 1277, 1361, 2123, 3791
	O,−1.685,−0.318,0.022
	O,0.954,1.087,0.172
	O,1.09,−1.160,−0.256
	H,1.612,−0.188,0.007
	H,1.205,−1.944,0.300
OSi(OH)2	Si,−0.070,−0.031,−0.076	8.466, 7.205, 3.892	313, 325, 358, 415, 478, 805, 865, 883, 988, 1296, 3851, 3855
	O,0.032,0.005,1.435
	O,1.154,0.011,−1.133
	O,−1.431,−0.119,−0.968
	H,0.929,−0.0217,−2.069
	H,−2.258,−0.151,−0.472
OSi(OH)2‐H2O complex	Si,−0.181,0.073,0.084	6.861, 1.632, 1.473	159, 185, 254, 296, 311, 316, 338, 377, 433, 618, 733, 805, 831, 902, 977, 1231, 1606, 3687, 3810, 3878, 3881
	O,0.223,1.331,−0.874
	H,0.388,2.158,−0.411
	O,−0.122,−1.218,−0.931
	H,0.112,−2.046,−0.499
	O,−0.019,−0.002,1.606
	O,−2.155,0.263,0.046
	H,−2.335,0.0820,0.986
	H,−2.593,−0.413,−0.493
TS from OSi(OH)2‐H2O to Si(OH)4	Si,−0.240,−0.003,0.168	4.602, 3.795, 3.409	983 i, 184, 264, 289, 317, 333, 375, 499, 524, 739, 770, 822, 826, 891, 1005, 1164, 1367, 2136, 3805, 3882, 3884
	O,−1.418,−0.597,−0.778
	H,−1.456,−1.554,−0.865
	O,−0.5442,1.598,0.234
	H,0,−0.227,2.047,1.024
	O,0.482,−0.791,1.316
	O,1.376,−0.216,−0.752
	H,0,1.506,−0.642,0.271
	H,1.911,0.571,−0.922
Si(OH)4	Si,−0.006,−0.127,0.163	4.114, 4.114, 3.654	189, 192, 281, 282, 285, 351, 365, 365, 399, 761, 845, 845, 847, 849, 932, 997, 997, 3886, 3886, 3387, 3890
	O,−0.259,0.016,1.775
	O,1.556,−0.245,−0.3159
	O,−0.731,−1.524,−0.291
	H,−0.260,−0.794,2.292
	O,−0.591,1.244,−0.516
	H,−1.374,1.631,−0.114
	H,0,2.048,0.577,−0.405
	H,−0.439,−1.919,−1.117

Importantly, both OSi(OH)2 and Si(OH)4 are stable with respect to reaction with O, H, and OH:

Thus, Si(OH)4 should be the main sink for meteor‐ablated silicon, as we demonstrate in section 3. Si(OH)4 should eventually polymerize with metallic species (e.g., FeOH, Mg(OH)2) to form MSPs. Although not the subject of this paper, we have shown that these polymerization reactions are exothermic. For example, the condensation reaction

Ion Chemistry

Metal atoms such as Fe and Mg become ionized in the lower thermosphere mostly through charge transfer with ambient NO+ and O2 + ions [Plane et al., 2015]. However, the concentration of Si must be negligible due to the fast reaction Si + O2 → SiO + O ((R1) in Table 1). The calculated adiabatic ionization energies of SiO and SiO2 are 11.5 and 12.5 eV, respectively. When compared with the calculated ionization energies of NO and O2, which are respectively 9.3 and 12.2 eV, it is clear that the only exothermic charge transfer reaction is

As shown in Figure 2, this appears to be the only route for ionizing neutral Si compounds. The resulting SiO+ must then be reduced to the observed Si+ (Figure 1) by reaction with atomic O:

Here we have assigned the rate coefficients for reactions (10) and (11) to their calculated Langevin capture rates ((R8) and (R9) in Table 1). It should be noted that Langevin theory provides an upper limit to the rate coefficient for an ion‐molecule reaction [Smith, 1980]. However, as is shown in the modeling section below, the unexpectedly high measured densities of Si+ relative to other metallic ions are already challenging to explain with both reactions (10) and (11) set to their Langevin capture rates, which implies that their rate coefficients must be very close to these capture rates.

Atmospheric Modeling

The global Si atmospheric model (termed here WACCM‐Si) is constructed from three components: (i) WACCM version 4, a “high‐top” coupled chemistry‐climate model with an upper boundary at 6.0 × 10−6 hPa (~140 km); (ii) a description of the neutral and ion‐molecule chemistry of Si (Figure 2 and Table 1); and (iii) the meteoroid input function (MIF), which specifies the temporal and geographical input of meteoric Si atoms into the MLT. The MIF was calculated using the astronomical characteristics of the sporadic meteor complex to estimate the global meteoric mass flux deposited at each altitude within the MLT, at all latitudes for each day of the year [Fentzke and Janches, 2008]. CABMOD was then used to compute the injection rates of individual metallic elements as a function of altitude [Vondrak et al., 2008]. Further details are given in our previous studies on Fe and Na [Feng et al., 2013; Marsh et al., 2013].

The Si MIF is illustrated in Figure 4. The peak Si injection rate of 0.018 atom cm−3 s−1 occurs around 97 km, and the annual mean Si ablation flux is 2.9 × 104 atom cm−2 s−1 (0.60 t d−1 over the Earth). It should be noted that this MIF predicts a Na flux of the required magnitude for modeling the mesospheric Na layer with WACCM [Marsh et al., 2013]. In fact, lidar measurements of the actual Na atom flux [Gardner et al., 2014; Huang et al., 2015] indicate that the Na flux may be underestimated by an order of magnitude, which would imply that WACCM requires additional vertical transport processes besides eddy diffusion in the MLT [Gardner and Liu, 2010]. From the point of view of the present study, where we are more concerned with the relative abundances of Si+, Fe+, and Mg+, this is less of an issue, although it does have implications for the transport of silicon species to high latitudes (see below).

Figure 4

(a) The seasonal variation of the silicon meteoric input flux (104 atom cm−2 s−1) as a function of latitude and month; (b) the annual average global Si injection profile.

WACCM‐Si was run from 2004 to 2013 (after a 10 year spin‐up), with dynamical fields nudged in the troposphere and stratosphere to the Goddard Earth Observing System 5 (GEOS5) meteorological data set, as described by Feng et al. [2013]. A nudging coefficient value (0.01) was used when assimilating the GEOS5 analysis into WACCM, so that 1% of the meteorological conditions were combined with WACCM fields below 60 km at every model dynamics time step. Above 60 km the model was free running. The WACCM‐Si output was sampled every 30 min for the rocket locations at 38°N and 68°N (see section 1).

Discussion

Figure 5 illustrates the vertical profiles of the major Si species predicted by WACCM‐Si at 38°N in August at local noon (11–13 LT). These profiles have been averaged from 7 years of data (2005–2011). A rocket‐borne profile of Si+ (flight 18.1006 [Herrmann et al., 1978]) for the same month, latitude, and local time is shown for comparison with the modeled Si+ profile. The layers of meteoric atomic ions tend to be highly perturbed by gravity waves and tides coupling with the Lorentz force from the Earth's magnetic field; these effects can be observed in the Ca+ layer, which (uniquely) has been observed with high spatial and temporal resolution by lidar [Gerding et al., 2000]. Since a rocket‐borne measurement provides a ~20 s snapshot, the agreement is surprisingly good. The sporadic E layer at 116 km is not predicted by WACCM (since these tend to be subscale, transient phenomena), but the underside of the Si+ layer is well captured by the model. Similarly good agreement is observed for flight 18.1008 at 37.8°N in January [Meister et al., 1978], and flight S37/P at 67.8°N in August [Kopp et al., 1984].

Figure 5

Modeled vertical profiles of the major neutral and ionic Si‐bearing species at 38°N, August and local noon. Profiles of Si+ (filled circles) and m/z 45, tentatively assigned to SiOH+ (grey squares), measured by rocket payload 18.1006 on 12 August 1976 [Herrmann et al., 1978].

Si+ is predicted to be the major silicon species only above 110 km, in contrast to metals such as Na and Fe where the ions dominate over the neutral atoms above ~95 km [Feng et al., 2013; Marsh et al., 2013]. There are two reasons for this: first, the relative difficulty of ionizing neutral silicon (which occurs only through charge transfer of SiO with O2 +, as discussed in section 2.2) and second, the reactions of Si+ and SiO+ with O3 ((R11) and (R15) in Table 1), which produce neutral SiO and SiO2, respectively. There are no counterparts to these reactions for the other meteoric metals. This explains the rapid decrease in Si+ below 97 km, illustrated in Figure 1a.

Figure 5 also shows that the only other significant ion is SiOH+, which forms through fast reactions of Si+ with H2O, SiO+ with H2 and H2O, and SiO2 + with H (Figure 2 and Table 1). Observed signals at m/z 45 [Goldberg and Aikin, 1972; Zbinden et al., 1975] have previously been assigned to SiOH+ [Kopp et al., 1995; Zbinden et al., 1975]. Although close to the background noise, the m/z 45 signal from payload 18.1006 shown in Figure 5 agrees very well at 105 km (45+/28+ = 0.004) with previous observations (45+/28+ = 0.003 [Goldberg and Aikin, 1972]) and is within a factor of 10 of the modeled [SiOH+] /[Si+] ratio at the Si+ peak. However, Figure 5 also shows an enhancement of the m/z 45 signal below 100 km, which is not consistent with the modeled profile of SiOH+. In view of the modeling results, it is unlikely that SiOH+, SiO2 +, and SiOH2 + contribute respectively to the m/z 45, 60, and 61 signals observed in other studies below 100 km and tentatively attributed to these species [Goldberg, 1975; Kopp et al., 1995].

Figures 1b and 1c compare the modeled [Si+]/[Fe+] and [Mg+]/[Fe+] ratio profiles, respectively, to the rocket averages (see section 1 for details of the five flights). Note that because the ratio from each of the five flights used to produce Figure 1 (see section 1) is averaged, the average ratios plotted in Figures 1b and 1c are different from the ratios of the averaged profiles in Figure 1a. The modeled ion ratio corresponding to each rocket flight is a 7 year average (2005–2011) value calculated within 2 h of the local time of the launch, and in the same month and location. The modeled Si+ is from this study, while the modeled Fe+ and Mg+ are from Feng et al. [2013] and Langowski et al. [2015], respectively.

The vertical dash‐dot‐dot lines in Figures 1b and 1c show the CI chondritic abundance ratio of Si and Mg to Fe [Asplund et al., 2009]. Above 105 km the [Si+]/[Fe+] ratio becomes superchondritic, which is indeed a common feature of all the rocket ion profiles reported to date (e.g., Aikin and Goldberg [1973], Goldberg [1975], Zbinden et al. [1975], and all the profiles used in this study from Kopp and coworkers). Comparison with the ratio calculated from WACCM shows that the model may underpredict the [Si+]/[Fe+] ratio, although the range of observed ratios overlaps the (much smaller) range of the modeled ratio. Inspection of Figure 5 shows that there is a substantial amount of SiO above 100 km, in excess over Si+ up to nearly 110 km. Although this is potentially an additional source of Si+, it should be remembered that the rate coefficients of reactions (9) and (11) are already set to their Langevin upper limits (Table 1) to maximize the production of Si+ from SiO via charge transfer with O2 +, and it is difficult to see another route to ionize SiO.

It is important to note that the Si MIF we have used here is not scaled with respect to the Na MIF. In contrast, when modeling the Fe and Mg layers using WACCM, the MIFs for Fe and Mg had to be divided by factors of 4 and 15, respectively, compared with Na [Feng et al., 2013; Langowski et al., 2015; Marsh et al., 2013]. In their modeling study, Fritzenwallner and Kopp [1998] also pointed out an apparent lack of differential ablation of Si with respect to Na. CABMOD treats meteoroids essentially as olivine particles (Fe‐Mg‐SiO4) which melt at 1850 K [Vondrak et al., 2008]. Si is predicted to ablate slightly less efficiently than Fe but has a slightly higher elemental abundance [Asplund et al., 2009]: the result is that for a range of MIFs with different mass and velocity distributions the Si/Fe ablation ratio ranges from 0.79 to 1.1 [Carrillo‐Sánchez et al., 2015]. In contrast, the implication of Figure 1b is that the Si/Fe ablation ratio is at least 4. One explanation for this may be that most of the cosmic dust particles entering the Earth's atmosphere is of cometary origin [Nesvorný et al., 2010], and there is some evidence that Fe is depleted by a factor of ~2 with respect to Si in particles returned from the Stardust mission [Flynn et al., 2006; Gainsforth et al., 2015]. However, that would not account fully for the seemingly efficient ablation of Si with respect to Fe.

WACCM‐Si correctly predicts the rapid decrease in the [Si+]/[Fe+] ratio below 100 km. Note that below 90 km, there is probably a significant contribution to the m/z 28 signal from H2CN+ [Kopp et al., 1995] (the 29Si/28Si ratio is also incorrect, but this may be due to a contribution from HCO+ at m/z 29 [Zbinden et al., 1975]). In contrast to Figure 1b, Figure 1c illustrates that WACCM captures very well the measured [Mg+]/[Fe+] ratio from 85 to 110 km, which is close to chondritic.

Figure 5 shows that SiO is predicted to be a major silicon species, occurring in a layer which peaks at ~97 km with a global average density of ~5000 cm−3. SiO is therefore the counterpart of the layers of metal atoms such as Na and Fe, which occur 5–12 km lower in the atmosphere. The SiO2 layer peaks around 96 km. SiO2 is a significant but never dominant species because it is formed via the oxidation of SiO by O3 and OH which are relatively slow ((R3) and (R4) in Table 1), compared with its reduction by H ((R5) in Table 1) or addition of H2O (R1). The global average lifetimes of Si+, SiO and SiO2 can be estimated by dividing the column abundances of their layers (Figure 5) by the meteoric ablation flux of Si (2.9 × 104 atom cm−2 s−1—see above). This yields individual lifetimes of 0.5, 1.9, and 0.6 days, respectively, or a total lifetime of these chemically labile Si species of 3.0 days.

Below 97 km the dominant sink is Si(OH)4. This species is then available to form MSPs by polymerizing with other metallic species. Viggiano et al. [1982] measured a negative ion with mass 76 and a density of a few 100 cm−3 below 80 km, which is about 2% of the predicted Si(OH)4 concentration below 80 km in Figure 5. Having argued that mass 76 ion was more likely SiO3 − than CO4 −, Viggiano et al. then speculated that SiO3 − could be produced by the rapid dissociative attachment of electrons to MSPs. Although recent rocket experiments [e.g., Plane et al. [2014]] have shown that electron attachment to MSPs produces heavy charged particles (>3000 amu), dissociative attachment to form SiO3 − could be a minor reaction channel.

Figure 6 illustrates the model diurnally averaged column abundances above 80 km of the major silicon species, Si+, SiO, and Si(OH)4, as a function of latitude and month (for the years 2004–2013). There are a number of striking features. First, Si+ and SiO exhibit a summertime maximum, particularly at high latitudes, where Si+ increases by a factor of ~10 and SiO by a factor of ~4, compared to midwinter. Second, the maxima in these species are not centered about the solstice but are shifted into late summer (by 2 months, in the case of SiO). Third, in contrast to Si+ and SiO, Si(OH)4 exhibits a wintertime maximum, again more pronounced at high latitudes where Si(OH)4 increases by a factor of ~5. Note that the maximum is in early winter, about 1 month before the solstice.

Figure 6

Diurnally averaged column abundances (109 cm−2) between 80 and 110 km of the major silicon species as a function of latitude and month, predicted by WACCM‐Si: (a) Si+, (b) SiO, and (c) Si(OH)4. Model output for the years 2004–2013 has been averaged.

The summertime maxima of Si+ and SiO are explained by the longer periods of sunlight causing increased ionization of the MLT, particularly at high latitudes. As shown in Figure 2, Si+ and SiO are closely coupled, principally through reactions (R8), (R9), (R10), (R11), and (R13) in Table 1. The shift in the peak toward autumn results from the Si MIF, which peaks in autumn and is a minimum in spring (Figure 4). The reason that a signature of the MIF appears in these species is that they are relatively short‐lived before conversion to the permanent sink Si(OH)4. This contrasts with the metals such as Na and Fe, where the metal atoms and ions persist much longer above 80 km before being permanently removed to MSPs; hence, long‐range transport tends to remove the signature of the MIF [Feng et al., 2013; Marsh et al., 2007]. Because Si(OH)4 forms rapidly via hydration of SiO2 (reactions (1) and (2)), and there is no sink apart from the slow formation of meteoric smoke, its global behavior reflects the meridional circulation in the upper mesosphere: high concentrations of Si(OH)4 arise from convergence of the meridional wind over the winter pole, particularly in early winter because of the autumnal maximum of the MIF in the relevant hemisphere, and low concentrations due to divergence of the meridional wind away from the summer pole.

Conclusions

In this study we present a new model for silicon chemistry in the MLT, based on a number of recent laboratory kinetic studies of Si, SiO, and Si+ [Eyet et al., 2010; Gómez Martin et al., 2009a; Gómez Martin et al., 2009b; Gómez Martín and Plane, 2011]. We have also used theoretical calculations to show that the likely fate of SiO2 in the MLT is a two‐step hydration to silicic acid. Si(OH)4 appears to be a very stable sink for silicon, whose likely fate is polymerization with metal oxides and hydroxides to form MSPs. The new model is therefore a significant advance on the previous work of Fritzenwallner and Kopp [1998].

Inclusion of this new silicon chemistry model in the 3‐D atmospheric model WACCM has revealed a number of interesting contrasts with other well‐known meteoric species such as Na and Fe. First, the counterpart of the Na and Fe layers is a layer of SiO, which is predicted to peak around 97 km. Second, the only route to ionize neutral silicon species appears to be charge transfer between SiO and O2 +, whereas metal atoms can also charge transfer with NO+ and undergo relatively fast photoionization [Plane et al., 2015]. The silicon oxide ions also undergo neutralizing reactions with O3 (Figure 2 and (R11) and (R15) in Table 1), which do not occur in the case of other meteoric species (e.g., Fe+ and FeO+). It is therefore particularly striking that Si+ has consistently been observed to be the major meteoric ion around 110 km. Since Fe+ and Mg+ should be the major forms of these two metals at 110 km [Feng et al., 2013; Langowski et al., 2015], this implies that the relative injection rate of Si from meteoric ablation, compared to metals such as Fe and Mg, is significantly larger than expected based on their relative chondritic abundances. This surprising observation is not yet fully explained. Finally, the global abundances of SiO and Si(OH)4 show clear evidence of the seasonal meteoric input function, which is seen to a much smaller extent in the case of other meteoric reservoir species such as FeOH [Feng et al., 2013], NaHCO3 [Marsh et al., 2013], and Mg(OH)2 [Langowski et al., 2015].