Model simulations of chemical effects of sprites in relation with
satellite observations

Abstract. Recently, measurements by the Superconducting Submillimeter-Wave Limb Emission Sounder (SMILES) satellite instrument have been presented which indicate an increase of mesospheric HO2 above sprite producing thunderstorms. These are the first direct observations of chemical sprite effects, and provide an opportunity to test our understanding of the chemical processes in sprites. In the present paper, results of numerical model simulations are presented. A plasma chemistry model in combination with a vertical transport module was used to simulate the impact of a streamer discharge in the altitude range 70–80 km, corresponding to one of the observed sprite events. Additionally, a horizontal transport and dispersion model was used to simulate advection and expansion of the sprite volumes. The model simulations predict a production of hydrogen radicals mainly due to reactions of proton hydrates formed after the electrical discharge. The net effect is a conversion of water molecules into H + OH. This leads to increasing HO2 concentrations a few hours after the electric breakdown. According to the model simulations, the HO2 enhancements above sprite producing thunderstorms observed by the SMILES instrument can not solely be attributed to the detected one sprite event for each thunderstorm. The main reason is that the estimated amount of HO2 released by a sprite is much smaller than the observed increase. Furthermore, the advection and dispersion simulations of the observed sprites reveal that in most cases only little overlap of the expanded sprite volumes and the field of view of the SMILES measurements is expected.


shows the key parameters of the measurements. In all three cases, the total enhancement of HO 2 is of the order of 10 25 molecules inside the field of view of the SMILES instrument. Note that the sprite locations lie outside the SMILES measurement volumes. The shortest horizontal distances between the SMILES lines of sight and the sprite bodies have been estimated to be about 10 km (events A and C) and 110 km (event B). There is a time lag of 1.5 to 4.4 hours between the 60 sprite detection and the SMILES measurement. Considering typical horizontal wind speeds in the upper mesosphere, Yamada et al. (2020) estimated advection distances of a few 100 km for the sprite air masses during the elapsed times between sprite occurrence and SMILES measurement. Data from the Worldwide Lightning Location Network (WWLLN) indicate strong lightning activity in the respective thunderstorm systems, and Yamada et al. (2020) pointed out that possibly additional sprites occurred which were not detected by ISUAL.

Model description
The main tool for our investigation is a one-dimensional atmospheric chemistry and transport model. It is used to simulate the undisturbed atmosphere before the occurrence of a sprite as well as the processes during and after the event. The model's altitude range is 40-120 km, and it's vertical resolution is 1 km. Table 2 shows the modelled species. The model's chemistry routines are based on a model which has previously been used to simulate short-term chemical effects of sprites (Winkler and 70 Notholt, 2014), and Blue Jet discharges (Winkler and Notholt, 2015). For a proper simulation of the atmospheric chemistry on longer time scales, the reaction scheme of this plasma chemistry model was merged with the one of an atmospheric chemistry model (Winkler et al., 2009) whose reaction rate coefficients were updated according to the latest JPL (Jet Propulsion Laboratory) recommendations (Burkholder et al., 2015). In the following, this model version is referred to as "Model_JPL". For some reactions also different rate coefficients have been considered. The reason is that there are reports on discrepancies between 75 modelled and observed concentrations of OH, HO 2 and O 3 in the mesosphere if JPL rate coefficients are used for all reactions (Siskind et al., 2013;Li et al., 2017). In particular, the JPL rate coefficient for the three body reaction where M denotes N 2 or O 2 , appears to be too small at temperatures of the upper mesosphere. According to Siskind et al. (2013), a better agreement between model and measurement is achieved if the rate coefficient expression proposed by Wong 80 and Davis (1974) is applied. Therefore, we have set up a model version "Model_WD" which uses the rate coefficient of Wong and Davis (1974) for reaction (1) while all other rate coefficients are as in Model_JPL. Li et al. (2017) have presented different sets of modified rate coefficients for reaction (1) as well as five other reactions of hydrogen and oxygen species. We have tested all these sets of rate coefficients in our model. In the following, we only consider the most promising model version which uses the 4th set of rate coefficients of Li et al. (2017). This model version is called "Model_Li4". Results of the model version 85 Model_JPL, Model_WD, and Model_Li4 will be compared in Sec. 4.
The effect of the enhanced electric fields occurring in sprite discharges is accounted for by reactions of energetic electrons with air molecules. The electron impact reaction rate coefficients are calculated by means of the Boltzmann solver BOLSIG+ (Hagelaar and Pitchford, 2005), for details see Winkler and Notholt (2014). For the present study, the electron impact reactions with H 2 O and H 2 shown in Tab. 3 have been added to the model.

90
The model has a prescribed background atmosphere of temperature, N 2 and O 2 altitude profiles. These profiles were derived from measurements of the SABER (Sounding of the Atmosphere using Broadband Emission Radiometry) instrument (Russell et al., 1999). The model uses daily mean day-time and night-time profiles calculated from SABER Level 2A, version 2.0, data for the geo-location of interest. At every sun rise or sun set event, the model's background atmosphere is updated.
The transport routines of the model calculate vertical transport due to molecular and eddy diffusion as well as advection. Details 95 on the transport model can be found in the Appendix A. Transport is simulated for almost all neutral ground state species of the model. Exceptions are N 2 and O 2 for which the prescribed altitude profiles are used. Transport is not calculated for ions and electronically excited species as their photochemical life-times are generally much smaller than the transport time constants.
The abundances of neutral ground state species at the lower model boundary (40 km) are prescribed using mixing ratios of a standard atmosphere (Brasseur and Solomon, 2005). At the upper model boundary (120 km), atomic oxygen and atomic 100 hydrogen are prescribed using SABER mixing ratios for an altitude of 105 km (The SABER O and H profiles extend only up to this altitude). This causes somewhat unrealistic conditions at the uppermost model levels, see Sec. 4. Furthermore, following Solomon et al. (1982), an influx of thermospheric NO is prescribed at the upper model boundary. For all other species a no-flux boundary condition is applied.

105
The one-dimensional model was used to simulate the atmosphere prior to sprite event B (Tab. 1). For this purpose, the model was initialised with trace gas concentrations from a standard atmosphere (Brasseur and Solomon, 2005) and then used to simulate a time period of almost ten years before the sprite event. The background atmosphere is made of zonal mean SABER profiles of the latitude stripe 0 • -13.5 • N (symmetric around the sprite latitude of 6.7 • N) for year 2009. For this spin-up run, no ionisation is included. This allows to use a chemical integration time step as large as one second. Transport is calculated once 110 every minute.
Three model versions with different vertical transport speeds have been tested, see Appendix A for details on the transport parameters. Figure 1 shows modelled mixing ratio profiles in comparison with SABER measurements as well as measurements by the MLS (Microwave Limb Sounder) satellite instrument (Waters et al., 2006). At high altitudes, above, say, 100 km, the modelled abundances of atomic oxygen and hydrogen are too small compared to SABER profiles. This is a result of using the SABER observations. Therefore, we decided to use this model version for the sprite simulations. The agreement between the model predictions and the measurements is not perfect but reasonable for a one-dimensional model compared to zonally averaged profiles.
The simulation results shown in Fig. 1 were obtained using JPL reaction rate coefficients (Model_JPL). The results of model  Model_WD agrees better with the SMILES HO 2 measurement, in particular if the vertical resolution of SMILES is taken into account (Fig. 2). On the other hand, Model_WD predicts too small OH number densities compared to the MLS measurements.

135
It is not possible to draw a conclusion here on what the best set of reaction rate coefficients is. We tend to favor Model_WD as it agrees better with the SMILES HO 2 data point than the other model versions. Note that the SMILES measurement corresponds to the actual geolocation and the solar zenith angle of the model profiles whereas the MLS data points are zonal night time averages. Furthermore, the MLS data had to be averaged over a wide latitudinal band (33 • S-40 • N) to reduce scatter or to obtain data points at all.

Sprite chemistry and vertical transport simulations
We have used the model version Model_WD with medium transport velocities for the sprite simulations. The sprite has been modelled as a streamer discharge at altitudes 70-80 km. This is the core region of the considered sprite (see plot (f) in the Supporting Information S1 of Yamada et al. (2020)). Following Gordillo-Vázquez and Luque (2010), a downward propagating streamer is modelled by an altitude dependent electric field time function which consists of two rectangular pulses, see Fig. 3.

145
The first pulse represents the strong electric fields at the streamer tip, the second pulse represents the weaker fields in the streamer channel (afterglow region). As in the study of Gordillo-Vázquez and , the electric field parameters are based on results of a kinetic streamer model (Luque and Ebert, 2010).
What follows here is a description of the model's parameters in terms of the reduced electric field strength E/N , where E is the electric field strength (V/cm) and N is the gas number density (cm −3 ). The reduced electric breakdown field strength is streamer head is reached. Based on the results of Luque and Ebert (2010), the streamer head peak electron density was assumed to scale with air density, and at 75 km altitude a value of 2 × 10 4 electrons per cm 3 is used.
Following Gordillo-Vázquez and , the reduced electric field of the second pulse is taken to be E k /N at all altitudes. Between the two pulses and up to a time of one second after the discharge, E/N has the sub-critical value 30 × 10 −17 Vcm 2 . The altitude dependent time-lag between the pulses is determined by the different propagation velocities of streamer head and streamer tail (Gordillo-Vázquez and . For the second pulse, a linearly decreasing duration 160 is assumed, ranging from 1.3 ms at 80 km to 0.3 ms at 70 km. As can be seen in Fig. 3, with this choice of parameters the modelled conductivities in the streamer channel are close to the value of 3 × 10 −7 (Ωs) −1 reported in the literature, see Gordillo-Vázquez and Luque (2010), and references therein.
The model was used to simulate a time period of 5.5 hours after the sprite event. It does not appear reasonable to simulate longer time periods with the one-dimensional model because eventually there is significant horizontal dispersion of the sprite 165 bodies, see Sec. 6. For the first two hours after the breakdown pulse, the full ion and excited species chemistry was simulated.
Then, the model switched into the less time-consuming mode without ion-chemistry (like in the spin-up run).
We begin our analysis of the chemical effects by an inspection of charged species. Figure 4 displays the simulated temporal evolution of the most important negative species at 80 km, and at 75 km. At lower altitudes, the general pattern is similar to the one at 75 km and is therefore not shown. The streamer tip electric field pulse leads to an increase of the electron density does not contribute significantly to the absolute electron detachment rates but plays a role for the hydrogen chemistry (see below). Subsequently, there is a formation of molecular ions, initiated mainly by e + O 3 → O − 2 + O at all altitudes. The resulting relative abundance of molecular ions is small at 80 km but larger at lower altitudes where eventually CO − 4 , CO − 3 , and Cl − become the most abundant ions. This is in overall agreement with previous model studies, e.g. Gordillo-Vázquez (2008); Sentman et al. (2008). 180 Figure 5 shows the simulated temporal evolution of the most important positive ions at 80 km, and at 75 km altitude. The primary ions resulting from electron impact ionisation of air molecules are N + 2 and O + 2 . The former undergoes rapid charge exchange (mainly) with O 2 , and after about one millisecond O + 2 has become the principal ion at all altitudes. This stays the same during the second electric field pulse. Eventually, there is a formation of O + 4 mainly through the three body reaction

185
The main loss process for O + 4 are reactions with water molecules: which produces a proton hydrate H + (H 2 O) 2 , and releases an OH radical. Larger proton hydrates can form via successive hydration: This is a well known mechanism in the D-region of the ionosphere (e.g., Reid (1977)), and was also predicted to take place in sprite discharges, e.g. Sentman et al. (2008); Evtushenko et al. (2013). According to our model simulations, proton hydrates  (3) and (7) are strongly pressure-dependent, and also because the abundance of water decreases with altitude ( Fig. 1).
Recombination of proton hydrates with free electrons release water molecules and atomic hydrogen: The same is true for recombination of proton hydrates with atomic or molecular anions. The net effect of the chain of reactions (3) -(8) is: As a result, there is a conversion of water molecules into two hydrogen radicals (H + OH).
Next we consider the impact on neutral hydrogen species. Figure 6 shows the decrease of water corresponding to the formation 205 of hydrogen-bearing positive ions and HO x at 75 km. The effect is similar at other altitudes 70-80 km. The small discontinuities of proton hydrates and HO x at a model time of two hours is due to the end of the ion chemical simulations (the total hydrogen amount is balanced, though). Figure 6 also displays the results of a model simulation of the undisturbed atmosphere, that is a model run without electric fields applied. Water, hydrogen radicals and several other species change in the no-sprite model simulation on time scales of hours (this is not a model drift but due to the fact that the night-time mesosphere is not 210 in a perfect chemical steady-state). Therefore, for a proper assessment of the sprite impact, the following analysis focuses on concentration differences between the sprite streamer simulation and the no-sprite simulation. Figure 7 shows the changes of the total amount of hydrogen atoms contained in those hydrogen species which are significantly affected by the sprite discharge at 75 km. During the first few seconds, there is a formation of positive hydrogen-bearing ions, and an increase of HO x at the expense of water molecules. After about ten seconds, the increase of HO x slows down, and water starts to recover.

215
The processes just analysed take place in the whole altitude range 70-80 km. At the highest altitudes, there is an additional process which affects the hydrogen chemistry. The already mentioned electron detachment process causes a conversion of H 2 into water molecules. This leads to a decrease of H 2 at 80 km compared to the no-sprite simulation, see Fig. 8. However, the production of HO x is still mainly due to hydration reactions of positive ions. The formation of HO x 220 molecules due to reactions of proton hydrates in the streamer at 80 km is smaller than it is at 75 km. The two main reasons for this are: (1) The total ionisation decreases with altitude (because the streamer tip peak electron density scales with air density); and (2)  The temporal evolution of the different HO x species at 75 km is resolved in Fig. 9. Initially, there is an increase of both OH and H concentrations due to the ion-chemical decomposition of water molecules while the concentration of HO 2 is decreased compared to the undisturbed atmosphere. The main reason for the latter are reactions of HO 2 with increased amounts of atomic oxygen produced in the sprite streamer: On longer time scales, the concentration of HO 2 increases. The most important production process at all altitudes is the three body reaction these parameters by considering the emissions in the first positive band of molecular nitrogen: there would be 2×10 21 excess molecules HO 2 . Even this value is small compared to the observed ∼10 25 molecules.

265
In conclusion, the estimated modelled enhancement of HO 2 due to a sprite is much smaller than the ∆HO 2 observed by the SMILES instrument. A possible reasons for these discrepancies could be missing chemical processes considered by the streamer model, inaccurate electric field parameters or reaction rate coefficients. All the results shown here were obtained with Model_WD. In order to test for the effect of changed rate coefficients, we have also performed sprite simulations using Model_JPL and Model_Li4. Generally, the results do not differ significantly from the Model_WD simulations. In particular, 270 the amount of HO 2 production is similar in all model versions. However, in case of Model_Li4, the formation of HO 2 is faster than in the other models. As a result, there is already an enhancement of HO 2 at a time of 1.5 hours after the electric breakdown pulse.
Another possible explanation for the differences between our model predictions and observations could be that there was an accumulation of HO 2 produced by a number of sprites. In the thunderstorm systems of interest, several hundreds of lightning 275 strikes occurred, and it is possible that multiple sprites developed (Yamada et al., 2020). Also, it is known that sprites tend occur The expansion of the sprite cross section is calculated by a Gaussian plume model approach, e.g. Karol et al. (1997). The radius of the plume corresponds to the standard deviation √ σ 2 of a Gaussian concentration distribution. If wind shear effects are neglected, the temporal change of the variance σ 2 is given by (Konopka, 1995): with K being the apparent horizontal diffusion coefficient. The formation time of a sprite is short compared to the time scales of atmospheric eddy diffusion. For such an instantaneous source, the diffusion coefficient is given by (Denison et al., 1994): where K ∞ is the atmospheric macroscale eddy diffusion coefficient, t is the age of the plume, and t L is the Lagrangian turbulence time scale. The latter is connected with K ∞ and the specific turbulent energy dissipation rate ε thought: Based on ranges of literature values of the turbulent parameters for the upper mesosphere (Ebel, 1980;Becker and Schmitz, 2002;Das et al., 2009;Selvaraj et al., 2014) we have considered two cases: (1) A slow diffusion scenario with K ∞ = 10 6 cm 2 s −1 , and ε = 0.01 Wkg −1 (2) A fast diffusion scenario with K ∞ = 2.5 × 10 7 cm 2 s −1 , and ε = 0.1 Wkg −1 315 10 https://doi.org/10.5194/acp-2020-1228 Preprint. Discussion started: 5 January 2021 c Author(s) 2021. CC BY 4.0 License.
The initial plume diameters were taken to be the horizontal widths of the sprites derived from the sprite observations (Tab. 1). Figure 12 shows results of the plume model simulations for both fast and slow diffusive sprite expansion. Only in case of sprite event C and for the fast expansion scenario, the SMILES field of view lies inside the expanded sprite body. For all other cases there is only little or no overlap of the SMILES measurement volume and the increased sprite volume. This is another indication that the measured HO 2 enhancements cannot (solely) be due to the three observed sprites. At all altitudes, the abundance of H 2 O is much larger than the produced amount of HO x . Therefore, water is not a limiting factor for the production of HO 2 . At all altitudes, the reaction H + O 2 + M → HO 2 + M is the most important process for the formation of HO 2 after the streamer discharge.
According to the model results, the HO 2 enhancements above sprite producing thunderstorms observed by the SMILES in-330 strument can not solely be attributed to the detected one sprite event for each thunderstorm system. The main reasons for this are: 1. The estimated amount of HO 2 released by a sprite is much smaller than the observed increase.

345
Appendix A: Transport modelling The transport part of the model calculates the change rate of the number density n i of species i according to the one-dimensional vertical diffusion and advection equation (Brasseur and Solomon, 2005, e.g.): with t being time, z altitude, D i the molecular or atomic diffusion coefficient of species i, K zz the vertical eddy diffusion 350 coefficient, α T the thermal diffusion factor, T the temperature in Kelvin, H i the individual scale height of species i, H the atmospheric scale height, and w the vertical wind speed. The diffusion coefficients D i (in cm 2 s −1 ) are given by (Banks and Kockarts, 1973): where M i and M are the molecular mass of species i and the mean molecular air mass (expressed in atomic mass units), 355 respectively, and n is the air number density (in units of cm −3 ). Following Smith and Marsh (2005), the thermal diffusion factor is taken to be α T = −0.38 for H and H 2 , and zero for all other species.
The free parameters in Eq. (A1) are the eddy diffusion coefficient K zz , and the vertical wind speed w. We have experimented with different altitude profiles of the eddy diffusion coefficient, and decided to use a profile parameterization proposed by 360 Shimazaki (1971): with standard coefficients S 1 = S 2 = 0.05 km −1 and S 3 = 0.07 km −1 . The parameter z 0 is the altitude at which the eddy diffusion is maximal, with K zz (z 0 ) = A. For all model simulations presented here, A = 10 6 cm 2 s −1 and z 0 = 105 km were used. The parameter B controls the eddy diffusion coefficient at lower altitudes. For the three cases "slow", "medium", and 365 "fast" vertical transport (Sec. 4), the following values have been used: B slow = 3 × 10 5 cm 2 s −1 , B medium = 5 × 10 5 cm 2 s −1 , and B f ast = 1 × 10 6 cm 2 s −1 .
There are one-dimensional model simulations of the middle atmosphere which do not consider vertical winds but only diffusive transport. We noted that the inclusion of advection due to winds significantly improves the model predictions compared to satellite measurements. In particular, the abundance of water in the middle to upper mesosphere increases and is in better to the diffusive transport. Therefore, the LIMA wind profile was multiplied by a scaling factor S < 1 to obtain the profile 375 Table 1. The three events analysed in Yamada et al. (2020). HW is the horizontal width of the sprites emissions, LT is the local time of the SMILES measurement, ∆t is the time difference between the sprite observation and the SMILES measurement, ∆r is the shortest distances between the field of view of the SMILES measurement and the sprite location, TH is the tangent height of the SMILES measurement, and ∆HO2 is the total enhancement along the line-of-sight of the SMILES measurement, *) Only a part of this sprite volume was observed by ISUAL (see Figure 1 in Yamada et al. (2020)). Table 2. Modelled species. The last row shows the molecules additionally included compared to the model of Winkler and Notholt (2015).
Negative species  Table 3. Electric field driven processes additionally included compared to the model of Winkler and Notholt (2014). The reaction rate coefficients in air depending on the reduced electric field strength were calculated with the Boltzmann solver BOLSIG+ (Hagelaar and Pitchford, 2005) using cross section data for H2O (Itikawa and Mason, 2005), and for H2 (Yoon et al., 2008).