Resolving the strange behavior of extraterrestrial potassium in the upper atmosphere

It has been known since the 1960s that the layers of Na and K atoms, which occur between 80 and 105 km in the Earth's atmosphere as a result of meteoric ablation, exhibit completely different seasonal behavior. In the extratropics Na varies annually, with a pronounced wintertime maximum and summertime minimum. However, K varies semiannually with a small summertime maximum and minima at the equinoxes. This contrasting behavior has never been satisfactorily explained. Here we use a combination of electronic structure and chemical kinetic rate theory to determine two key differences in the chemistries of K and Na. First, the neutralization of K+ ions is only favored at low temperatures during summer. Second, cycling between K and its major neutral reservoir KHCO3 is essentially temperature independent. A whole atmosphere model incorporating this new chemistry, together with a meteor input function, now correctly predicts the seasonal behavior of the K layer.


Introduction
The ablation of meteoroids in the upper atmosphere leads to the formation of distinct layers of metal atoms around 90 km. These layers, which are global in extent, provide a sensitive marker for the dynamics and chemistry of the atmosphere at the edge of space [Plane, 2003]. The presence of K atoms in the mesosphere and lower thermosphere (MLT) was first established from photometric observation of the twilight emission at 766.5 and 769.9 nm [Sullivan and Hunten, 1964]. Subsequent work using the resonance lidar technique [Megie et al., 1978;Eska et al., 1999] showed that the K column density varies semiannually, with maxima in summer and winter and minima at the equinoxes. The layer thus exhibits little overall seasonal variation at midlatitudes. In contrast, the Na layer has a pronounced annual variation, being around 3 times larger in winter than summer . As a result, the Na/K abundance ratio at midlatitudes decreases from~50 in winter to only~10 during summer. The explanation for this very surprising difference in behavior between two Group 1 (alkali) metals, whose chemistry should be very similar, is a problem of 50 years' standing.
Initially, attention focused on whether the Na/K ratio indicated a seasonal change in the source of the metals in the MLT. The Na/K elemental ratio in meteoroids is around 15, whereas the seawater ratio is 47, suggesting that the vertical transport of sea-salt particles during storms in the polar night could be responsible for the wintertime enhancement of Na over K, which is greatest at high latitudes [Sullivan and Hunten, 1964;Megie et al., 1978]. However, in addition to the challenge of transporting material upward within the winter polar vortex to the MLT where the mean wind is downward, Jegou et al. [1980] showed that the mesospheric Na/Li ratio remains constant or may even decrease during winter, which appears to rule out a sea-salt source. An alternative explanation proposed that the neutral metal atom chemistry in the MLT is controlled by the dynamics of hydrated metal ion clusters [Jegou et al., 1985], but the assumed temperature-dependent differences in the ion cluster chemistries required for this model have not subsequently been validated by experiment or theory. More recent modeling attempts [Eska et al., 1999;Delgado et al., 2006] have been unable to identify key differences in the chemistry of K and Na. K + ions have been observed by rocket-borne mass spectrometry in a permanent layer above 90 km [Kopp, 1997]. These ions form continuously from K in the MLT via photoionization and charge transfer with ambient ions [Eska et al., 1999]: ©2014. The Authors.

RESEARCH LETTER
Neutralization of K + occurs by forming a cluster ion with a ligand X, where X = N 2 , O 2 , O (in order of decreasing concentration in the MLT). These clusters can then undergo ligand switching with another atmospheric component Y (CO 2 or H 2 O) to form a more stable cluster, before undergoing dissociative electron recombination: where M is a third body (most likely N 2 or O 2 ). Plane et al. [2006] carried out high-level quantum chemistry calculations on the K + · X and · X · K + · Y cluster ions and then used Rice-Ramsperger-Kassel-Markus (RRKM) theory to estimate the rate coefficients for reactions (4) and (5). The important result is that because K + is a relatively large singly charged ion (cf. Na + ), it forms very weakly bound clusters with the major species (N 2 , O 2 , and O) in the MLT, with binding energies less than 20 kJ mol À1 . The sequence of reactions (4)-(6) therefore only becomes significant at the very low temperatures characteristic of the summertime MLT. Figure 1 illustrates schematically the important pathways that couple K + and K in the MLT.
Although the ion-molecule chemistry described above should produce an enhancement of K during summer, neutral chemistry would be expected, by analogy with Na [Plane, 2004], to remove K through formation of stable reservoir species. In the case of Na, the important reservoir is sodium bicarbonate (NaHCO 3 ), which is only converted back to atomic Na by photolysis or reaction with H. The reaction NaHCO 3 + H → Na + H 2 CO 3 has an activation energy of about 10 kJ mol À1 [Cox et al., 2001]. Thus, during summer the reaction becomes very slow, causing NaHCO 3 to build up and contributing to the observed minimum in the Na column density [Plane, 2004].
In this paper we use a combination of theoretical techniques to explore the analogous chemistry of potassium. The resulting rate coefficients for the neutral reactions, together with those from our previous study of K ion-molecule reactions [Plane et al., 2006], are then incorporated into a whole atmosphere chemistry-climate model.

Theoretical Calculations on Neutral K Chemistry
In order to calculate rate coefficients and photolysis rates, a variety of electronic structure calculations was performed using the Gaussian 09 suite of programs [Frisch et al., 2009]. Initially, the hybrid density functional/ Hartree-Fock B3LYP method was employed with the 6-311 + G(2d,p) triple zeta basis set, which has both polarization and diffuse functions added to the atoms. Excited state calculations were performed using the time-dependent TD/B3LYP method [Bauernschmitt and Ahlrichs, 1996]. Very accurate energies of the stationary points (i.e., local minima and transition states) on the relevant potential energy surfaces (PESs) were calculated at the complete basis set (CBS-QB3) level of theory [Montgomery et al., 1999]. In the case of the KO radical, the CBS-QB3 method failed due to excessive mixing of the frozen core and valence orbitals in the coupled cluster step of the calculation. Accurate reaction enthalpies for reactions involving KO were therefore carried out at the B3LYP level using the augmented version of the relativistic small-core effective core potential ECP10MDF basis set of the Stuttgart/Cologne group [Lim et al., 2005]. As shown in Figure 1, K atoms are oxidized by O 3 to yield KO, which a second O 3 can then oxidize to KO 2 [Plane, 2002]. KO 2 also forms directly from the recombination of K and O 2 [Plane et al., 1990]. However, in the MLT these oxides exist in the presence of O, H, H 2 , H 2 O, and CO 2 . Calculations (see above) show that the following reactions are all exothermic: Reaction (9) has a second channel forming K + H 2 O which is exothermic by 228 kJ mol À1 . However, in the case of the analogous NaO + H 2 reaction, formation of Na is the minor channel [Ager and Howard, 1987] and we assume the same for K. Reaction (11) could also form KOH + O, which is exothermic by 135 kJ mol À1 but involves breaking the strong O-O bond in the superoxide and is therefore probably slow.
The potential energy surfaces of reactions (7)-(12) do not contain significant energy barriers, and so by analogy with the corresponding reactions of NaO and NaO 2 [Ager and Howard, 1987;Helmer and Plane, 1993;Griffin et al., 2001], these reactions are likely to be fast with little temperature dependence. We therefore assign them a typical collision theory rate coefficient of 2 × 10 À10 exp(À120/T) cm 3 molecule À1 s À1 .
KOH can also recombine with CO 2 to form the bicarbonate: Figure 2a illustrates the potential curve for reaction (13), together with the geometries of the stationary points (atom coordinates are listed in Table S1 in the supporting information). The KOH and CO 2 initially recombine to form the K(OH)CO 2 adduct. This adduct then rearranges to the more stable KHCO 3 either via migration of the H or the K atom, the latter having a much lower energy barrier with respect to the energy of the reactants. The rate coefficient for reaction (13) was calculated using a multiwell energy-grained master equation solved with the program Master Equation Solver for Multienergy Well Reactions [Robertson et al., 2008], where the microcanonical rate coefficients linking the reactants, intermediates, and products were calculated by RRKM theory using the molecular parameters in Table S1 (supporting information). The internal energies of the intermediates on the PES were divided into a contiguous set of grains (width 100 cm À1 ), each containing a bundle of rovibrational states. For dissociation to products or reactants, microcanonical rate coefficients were determined using inverse Laplace transformation to link them directly to the dipole-induced dipole capture rate coefficient [Georgievskii and Klippenstein, 2005], k capture (KOH + CO 2 ) = 5.8 × 10 À10 (T/300 K) 1/6 cm 3 molecule À1 s À1 . The probability of collisional transfer between grains was estimated using the exponential down model (average energy for downward transitions was set to < ΔE > down = 300 cm À1 , assumed temperature independent). At the low-pressure limit appropriate for the MLT region, k 13 (140-220 K) = 7.1 × 10 À28 (200 K/T ) 4.21 cm 6 molecule À2 s À1 .
KHCO 3 could potentially be destroyed by reaction with H or O: These reactions are either exothermic or only slightly endothermic and so were investigated in more detail. Figure 2b shows the potential energy curves and transition state geometries for the three channels of reactions (14a)-(14c). Channels (14a)-(14c) have barriers of 32, 61, and 71 kJ mol À1 , respectively. Channels (14b) and (14c) can therefore be ruled out at MLT temperatures. The rate coefficient for reaction (14a) was estimated to be k 14a (140-220 K) = 4.5 × 10 À11 exp(À3590/T) cm 3 molecule À1 s À1 by using transition state theory with a correction for tunneling of H through the reaction barrier, assumed to be a symmetrical Eckart barrier as in our previous study of the analogous reaction of NaHCO 3 [Cox et al., 2001] (see Table S1 in the supporting information for the relevant molecular parameters). At the various levels of theory used in this study, the energy barrier for KHCO 3 + H is consistently between 7.1 and 9.3 kJ mol À1 higher than that for NaHCO 3 + H. Thus, even allowing for some uncertainty in the calculated absolute value of the barrier for reaction (14a), the reaction should be at least 2 orders of magnitude slower at 200 K than its Na analogue. Reaction (15) has a significant barrier of 38 kJ mol À1 and will therefore be too slow in the MLT for further consideration.
The photolysis cross sections for KHCO 3 , KO 2 , and KOH as a function of wavelength were calculated at the TD/ B3LYP/6-311 + g(2d,p) level of theory [Frisch et al., 2009] and then placed on an absolute scale using our previously measured cross sections for the Na analogues [Self and Plane, 2002]. The photolysis rates of these three potential molecular reservoirs were estimated to be 1.2 × 10 À4 , 2.2 × 10 À3 , and 0.027 s À1 , respectively. Figure S1 (supporting information) illustrates the cross-section spectra for KHCO 3 and NaHCO 3 , which have thresholds around 280 nm.
Finally, the dimerization rate coefficient of KHCO 3 was estimated from the dipole-dipole capture rate [Georgievskii and Klippenstein, 2005] to be 1.0 × 10 À9 (T/200 K) À0.29 cm 3 molecule À1 s À1 , using a calculated value for the dipole moment of KHCO 3 of 8.2 Debye (at the B3LYP/6-311 + g(2d,p) level). However, KHCO 3 can also polymerize with other meteoric constituent molecules (e.g., NaHCO 3 , FeOH, and Mg(OH) 2 ), and the dimerization rate coefficient needs to be increased to account for this. The meteoritic elemental ratios of Na, Fe, Mg, and Si relative to K are 15.5, 248, 279, and 264, respectively [Asplund et al., 2009]. Thus, the concentration of compounds containing these other elements could potentially be~800 times higher than the concentration of KHCO 3 . However, recent modeling studies [Vondrak et al., 2008;Feng et al., 2013;Langowski et al., 2014] indicate that Fe, Mg, and Si ablate with much lower efficiencies than K (and Na), so the polymerization rate coefficient should be estimated by increasing the dimerization rate coefficient by significantly less than 800. A factor of 270 (i.e., a phenomenological polymerization rate coefficient of 2.7 × 10 À7 cm 3 molecule À1 s À1 ) was found to reproduce the observed underside of the K layer satisfactorily, as shown below. Table 1 contains a complete list of rate coefficients for the neutral and ion-molecule reactions depicted in Figure 1. In order to model the K layer, we adopt the same approach that was used recently to model the global Na, Fe, and Mg layers Marsh et al., 2013;Langowski et al., 2014]. Three components are combined to construct the K model (termed here WACCM-K): (i) the Whole Atmosphere Climate Community Model (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 K ( Figure 1 and Table 1) WACCM-K was run from 2004 to 2012 (after a 9 year spin-up), with dynamical fields nudged in the troposphere and stratosphere to the Goddard Earth Observing System-5 (GEOS-5) meteorological data set, as described by Feng et al. [2013]. A nudging coefficient value (0.01) was used when assimilating the GEOS-5 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 MIF is calculated from current knowledge of the astronomical characteristics of the sporadic meteor complex to estimate the
The WACCM-K output was compared with lidar measurements of the K layer made over 3 years (2004)(2005)(2006) at Kühlungsborn, Germany (54°N, 12°E) (see Eska et al. [1999] for a description of the technique). The K layer at this station is typical of other midlatitude locations [Sullivan and Hunten, 1964;Megie et al., 1978]. To minimize the impact of solar cycle effects on the K layer, we focus here on the seasonal variation of the K layer averaged over the 3 years. Figure 3a shows the measured seasonal variation of the K density as a function of height, and Figure 3b shows the modeled variation. The first noteworthy point is that the semiannual variation of the K layer, with a midsummer maximum and minima at the equinoxes, is very well captured by the model. A second point is that the absolute K density agrees well with the observed density, using a MIF for K (average 120 atom cm À2 s À1 ) which is in the chondritic ratio (K/Na = 6.5% [Asplund et al., 2009]) to the Na MIF used in our recent study of the Na layer . This is in accord with a chemical ablation model [Vondrak et al., 2008] which predicts that Na and K should ablate with similarly high efficiencies compared to other meteoric metals. A third point is that the peak height and width of the layer are satisfactorily captured by the model (see Figure S2 in the supporting information which shows MLT profiles of all the K-containing species in WACCM-K). Figure 3c shows that the modeled semiannual seasonal variation of the K column abundance is also in good agreement with the lidar measurements. The contrasting seasonal variation of the Na layer at the same location  is included for comparison.

Discussion and Conclusions
Why does K exhibit a summertime peak and equinoctial minima, in contrast to Na? There are two factors involved. First, the K + ion is more difficult to neutralize because it can only form weakly bound cluster ions at very low temperatures. Thus, the low temperatures of the summertime MLT shift the balance from K + to K, accounting for the summertime maximum. The second factor concerns the neutral chemistry. None of the reactions involving KO, KO 2 , and KOH has a significant temperature dependence. In contrast, the activation ) as a function of altitude and month, measured by lidar. (b) K density as a function of altitude and month from WACCM-K. (c) The observed and modeled K column abundance, compared with the Na column abundance (since Na lidar measurements are not available from Kühlungsborn for the same period, validated model output  for that location averaged over 2004-2006 is shown). energy for the reaction KHCO 3 + H is so large that this reaction never plays a significant role, even during the warmest period in winter. The only route from KHCO 3 back to K is via photolysis, and this is essentially temperature independent. For Na, the reaction NaHCO 3 + H has a smaller activation energy [Cox et al., 2001], and so during winter this reaction becomes fast enough to recycle the reservoir to atomic Na, leading to the pronounced wintertime maximum. Because there is less winter-summer variation in the K layer, the semiannual variations become pronounced. As shown in Figure S3 (supporting information), atomic O exhibits minima at the equinoxes, whereas H 2 and CO 2 exhibit maxima. Inspection of Figure 1 shows that simultaneously decreasing O and increasing H 2 and CO 2 will drive K to KHCO 3 , and this explains the minima in the K layer seen during the equinoxes.
It should be pointed out that the MIF used here would correspond to a global input rate of 4.6 t d À1 of ablated cosmic dust . This is at the lower end of current estimates of this quantity, which vary from~3 to 300 t d À1 [Plane, 2012]. Although this estimate is calculated from an astronomical model constrained to measurements of meteor head echoes using a number of high-performance large aperture radars [Fentzke and Janches, 2008], it is possible that a population of undetected meteors could make a significant contribution to the MIF. To test this possibility, we ran a sensitivity experiment where the MIF was increased by a factor of 10 in WACCM-K. The K column density increased by a factor of 8.6 (there is a small nonlinearity because of the second-order polymerization reaction R13 in Table 1). The coefficient of R13 would therefore need to be increased if the MIF was larger than 4.6 t d À1 . However, the layer height, the layer shape, and the semiannual seasonal behavior of the K layer-which is the focus of this study-were unchanged.
In conclusion, the surprisingly different behavior of the K layer now appears to be explained. Nevertheless, it would be desirable if the theoretical rate calculations presented here for reactions (13) and (14a)-(14c), and the assumed rate coefficients for reactions (8)-(12), could be tested experimentally in the future.