D region ion-neutral coupled chemistry within a whole atmosphere chemistry-climate model

This study presents a new ion-neutral chemical model coupled into the Whole Atmosphere Community Climate Model (WACCM). The ionospheric D region (altitudes ~50 90 km) chemistry is based on the Sodankylä Ion and Neutral Chemistry (SIC) model, a 1-dimensional model containing 306 ion-neutral and ion-recombination reactions of neutral 15 species, positive and negative ions, and electrons. The SIC mechanism was reduced using the Simulation Error Minimization Connectivity Method (SEM-CM) to produce a reaction scheme of 181 ion-molecule reactions. This scheme describes the concentration profiles at altitudes between 20 km and 120 km of a set of major neutral species (HNO3, O3, H2O2, NO, NO2, HO2, OH, N2O5) and ions (O2 + , O4 + , NO + , NO + (H2O), O2 + (H2O), H + (H2O), H + (H2O)2, H + (H2O)3, H + (H2O)4, O3 , NO2 , O , O2, OH , O2 (H2O), O2 (H2O)2, O4 , CO3 , CO3 (H2O), CO4 , HCO3 , NO2 , NO3 , NO3 (H2O), NO3 (H2O)2, NO3 (HNO3), NO3 20 (HNO3)2, Cl , ClO ), which agree with the full SIC mechanism within a 5% tolerance. Four 3D model simulations were then performed, using the impact of the January 2005 Solar Proton Event (SPE) on D region HOx and NOx chemistry as a test case of four different model versions: the standard WACCM (no negative ions and a very limited set of positive ions); WACCM-SIC (standard WACCM with the full SIC chemistry of positive and negative ions); WACCM-D (standard WACCM with a heuristic reduction of the SIC chemistry, recently used to examine HNO3 formation following an SPE); and 25 WACCM-rSIC (standard WACCM with a reduction of SIC chemistry using the SEM-CM Method). Standard WACCM misses the HNO3 enhancement during the SPE, while the full and reduced model versions predict significant NOx, HOx and HNO3 enhancements in the mesosphere during solar proton events. The SEM-CM reduction also identifies the important ionmolecule reactions that affect the partitioning of odd nitrogen (NOx), odd hydrogen (HOx), and O3 in the stratosphere and mesosphere. 30 Geosci. Model Dev. Discuss., doi:10.5194/gmd-2016-57, 2016 Manuscript under review for journal Geosci. Model Dev. Published: 20 April 2016 c © Author(s) 2016. CC-BY 3.0 License.

mechanism was reduced using the Simulation Error Minimization Connectivity Method (SEM-CM) to produce a reaction scheme of 181 ion-molecule reactions of 181 ion-molecule reactions of 27 positive and 18 negative ions.This scheme describes the concentration profiles at altitudes between 20 km and 120 km of a set of major neutral species (HNO 3  ), NO 3 -(HNO 3 ) 2 , Cl -, ClO -), which agree with the full SIC mechanism within a 5% tolerance.Four 3D model simulations were then performed, using the impact of the January 2005 Solar Proton Event (SPE) on D region HO x and NO x chemistry as a test case of four different model versions: the standard WACCM (no negative ions and a very limited set of positive ions); WACCM-SIC (standard WACCM with the full SIC chemistry of positive and negative ions); WACCM-D (standard WACCM with a heuristic reduction of the SIC chemistry, recently used to examine HNO 3 formation following an SPE); and WACCM-rSIC (standard WACCM with a reduction of SIC chemistry using the SEM-CM Method).Standard WACCM misses the HNO 3 enhancement during the SPE, while the full and reduced model versions predict significant NO x , HO x and HNO 3 enhancements in the mesosphere during solar proton events.The SEM-CM reduction also identifies the important ion-molecule reactions that affect the partitioning of odd nitrogen (NO x ), odd hydrogen (HO x ), and O 3 in the stratosphere and mesosphere.

Introduction
Energetic charged particles that impact on the Earth's atmosphere come from several different sources: in the case of protons (and some heavier ions), directly from the Sun during Solar Proton Events (SPEs) and from outside the Solar System in the form of high-energy Galactic cosmic rays; and in the case of electrons, from the radiation belts around the Earth during geomagnetic storms and sub-storms (Du et al., 2008;Kurt et al., 2004;Shea and Smart, 1992;Sinnhuber et al., 2012).
Energetic particle precipitation (EPP) can lead to significant perturbations of the chemical composition from the lower thermosphere all the way down to the stratosphere (Jackman et al., 2014;Verronen et al., 2011).Protons with energy above 1 MeV can penetrate down to the mesosphere and the upper stratosphere, particularly at high geomagnetic latitudes.
Energetic Particle Precipitation (EPP) causes ion-pair formation, and the subsequent neutralization produces odd nitrogen (NO x = N + NO + NO 2 ) and odd hydrogen (HO x = H + OH + HO 2 ) species.The NO x and HO x species destroy mesospheric ozone via catalytic cycles (Crutzen, 1970;McElroy et al., 1992), which can have a significant impact on the radiative balance of the middle atmosphere and hence on climate (Sinnhuber et al., 2012).
Current whole atmosphere chemistry climate models such as the Hamburg Model of the Neutral and Ionized Atmosphere (HAMMONIA) (Schmidt et al., 2006) and the Whole Atmosphere Community Climate Model (WACCM) (Garcia et al., 2007;Jackman et al., 2011;Marsh et al., 2007;Marsh et al., 2013) have simple lower E region plasma chemistry, essentially NO + and O 2 + ions balanced by free electrons, but do not describe the much greater complexity of the D region where clusters and negative ions dominate (Brasseur and Solomon, 2005;Sinnhuber et al., 2012;Winkler et al., 2008).Therefore, in order to model the impacts of EPP in the atmosphere below 90 km, it is essential to have a detailed treatment of D region ion chemistry.
One of The the leading kinetic models of D region chemistry is the Sodankylä Ion and Neutral Chemistry (SIC) model, developed jointly at the Sodankylä Geophysical Observatory and the Finnish Meteorological Institute (Turunen et al., 1996;Verronen et al., 2011).It contains 306 307 ion-molecule and also 2254 ion-ion recombination reactions.The model is a 1D model with simple vertical transport (molecular and eddy diffusion) to balance the computational cost.Putting such a large number of additional reactions into a 3D global model could be computationally expensive because of the large number of stiff partial differential equations that need to be solved.In view of this, several of the authors of the present paper carried out a heuristic reduction of the full SIC chemistry scheme -i.e., based on chemical knowledge and intuition -in order to produce a D region ion scheme of 196 ion-neutral reactions, which was able to account for the substantial enhancements of HNO 3 that have been observed between ~50 and 80 km after major SPEs (Andersson et al., 2016).
Encouraged by the development of Whole Atmosphere Community Climate Model with D region chemistry (WACCM-D), we have now carried out a systematic mechanism reduction of the SIC chemistry, which is the subject of this paper.The objective of the reduction was to model a set of important species during a highly perturbed period -the intense SPE of late October 2003 -at an acceptable level of accuracy, but with reduced computational time.The reason for choosing a highly perturbed period is that with this it was possible to reveal every important chemical detail, while using a weaker event would have caused some important reactions to be missed.The full and reduced SIC model chemistries were then tested in WACCM, where the medium SPE of 15-17 January 2005 was used as a test case of how well the reduced scheme captures the substantial atmospheric perturbation to NO x , HO x and HNO 3.

Methodology
The SIC model (Verronen et al., 2011) contains 306 307 ion-neutral reactions in total of all the important and necessary species.These are listed in the Supplementary Material (SM), where the reactions in the grey shaded rows are the subset selected systematically in this study for the reduced model, and the reactions shown in red bold type are in WACCM-D (Andersson et al., 2016).It should be noted that rate coefficients for some of the ion-molecule reactions have been updated using a recent review (Pavlov, 2014) -these changes are documented in the SM.

Description of the 1D SIC model
In the initial step, primary proton collisions with molecular nitrogen or oxygen cause dissociative ionization and produce secondary electrons which also form atomic N or O from N 2 and O 2 .A detailed description of the SIC model is given in (Verronen et al., 2005), where the ionization scheme is also described.It should be noted that the SIC model assumes that the electron temperature is the same as that of the neutral atmosphere, which is a reasonable simplification in the D region (Roble, 1976).The model solves time-dependent concentrations between 20 km and 150 km with 1 km vertical resolution.
The concentrations are controlled by solar radiation, particle precipitation, chemical reactions and vertical transport.Daily solar irradiance data are taken from the SOLAR2000 model (Tobiska and Bouwer, 2006).The integrated proton fluxes measured by the Geostationary Operational Environmental Satellite (GOES-11) satellite are converted to energy resolved flux spectra using the exponential rigidity relation (Freier and Webber, 1963).Ionization rates are calculated from the spectra based on the proton energy-range measurements in standard air (Bethe and Ashkin, 1953) as described in (Verronen et al., 2005), assuming that 35 eV of energy is required to produce one ion pair (Porter et al., 1976) , 15, 20, 25, 30, 45, 60 and 100 km are taken from (Shimazaki, 1984).

Mechanism reduction
The updated SIC mechanism was reduced using the Simulation Error Minimization Connectivity Method (SEM-CM) (Nagy and Turanyi, 2009).SEM-CM is based on the Connectivity method (Turanyi, 1990;Turanyi et al., 1989), which identifies necessary species based on their kinetic connectivity to the designated important species.The important species are those whose concentration should be accurately reproduced by the reduced model.Further necessary species are those whose inclusion in the reduced mechanism is also needed to achieve this aim.
As a result, the reduced mechanism contains only the reactions of the important and the necessary speciesthat is, the nonredundant reactions.The set of redundant species take part only in the redundant subset of reactions.In addition, in the mechanism reduction process the redundant species are identified.In the following, redundant reactions refer to the set of reactions involving redundant species.
Kinetic connectivity is defined via the kinetic differential equation which describes the change of species concentrations (vector f = (f 1 ,...,f N ) for N species) due to thermal ( number of reactions) and photochemical ( number of reactions) transformations: Vector function f defines the rate of concentration change of all species as a function of the concentration of all species (c=(c 1 ,...,c N )) and the thermal and photochemical rate parameters (vectors k=(k 1 ,...,k  ) and J=(J 1 ,...,J  ) respectively), which depend on temperature (T), pressure (p) and solar actinic flux (function F(), where  is wavelength).The dimension of J is 23, which equals the number of photodissociation ( 16) and photoionization ( 7) reactions (these are listed in Table 2), while the dimension of k is the total number of neutral-neutral and ion-neutral thermal reactions (306307).Note that the photoionization and photodissociation reactions are not listed in the SM as no fixed rate parameters are defined for them (they depend on the solar activity).The kinetic connectivity B j of species j to the set of i selected species is measured by the sum of the squared elements of the log-normalized Jacobian matrix ( ij J = log f i / log c j , where matrix J is an NN matrix) of the kinetic differential equation.
Species with relatively large B j values are closely linked to the set of selected species (which initially are the important species), and their presence is necessary in the mechanism for the formation and consumption of the selected species.The Jacobian matrix is determined at several time points from a simulation with the full model, and then stored for the iteration.
The reduced mechanism is constructed in an iterative procedure by adding additional species to the set of selected species in order to achieve the stipulated accuracy with which the concentrations of the important species need to be reproduced.Note that the initial set of necessary species is the set of important species for which the concentrations are to be reproduced within the given accuracy.
However, adding a single species to the important ones will not necessarily require the inclusion of additional reactions, thus species should be added in so-called complementary sets.A complementary set contains all species from a reaction which has not yet been selected.The kinetic connectivity of the k-th complementary set (C k ) with n k species to the selected species is defined as the average of the B values: At each time of the investigation, the complementary sets are ranked according to their connectivity and several reduced mechanisms are formed by adding the strongly connected sets to the selected species and including their reactions in the reduced scheme.The SEM-CM mechanism at this step investigates whether each species in the reduced mechanism is either an initially present species or there is a chemical pathway of formation from the species initially present.Additional complementary sets are added until this condition is fulfilled (Nagy and Turanyi, 2009).The mechanism reduction is carried out in a gradual species inclusion procedure ("building"), and in each step of this, several candidate extended sets of species and corresponding reduced mechanisms are generated.The candidate reduced mechanisms are simulated and stored with their simulation errors in a database sorted by their numbers of species.In the database the mechanisms with one more species and the one with smallest simulation error is selected.If the mechanism fulfils the required accuracy criterion, the reduction procedure is complete.If this is not the case, the SEM-CM building step is repeated with these "selected" species of smallest error until the required accuracy is met.In essence, the construction of the reduced mechanism according to the SEM-CM procedure is governed by the steepest decrease of simulation error.
The aim of the mechanism reduction procedure is to reduce the maximum error to below a defined threshold.To achieve this, simultaneous minimization of a maximum ( MAX ) and a root-mean-square global error ( RMS ) by two SEM-CM threads was found to be very efficient.These errors were defined as: This mixed error behaves as a relative error close to the maximum concentration of species j in scenario i ( full MAX , j i c ) obtained with the full mechanism, whereas at much lower concentrations it damps large relative deviations and behaves as a scaled absolute error.
For the SIC mechanism reduction, the maximum of a large SPE was selected in order to be able to ensure that the reduced scheme could satisfactorily model large, short-lived perturbations to ions and neutrals.This was the SPE in October 2003, the "Halloween Storm" (Burlaga et al., 2005) that peaked at 0615 UT on 29 th October 2003 with a proton flux of 29,500 (pfu>10MeV) (pfu: particle flux unit = particle cm -2 ster -1 s -1 ).Four specific altitudes, 60, 70, 80 and 90 km were selected to represent the D region.The selection of important chemical species contained those that are the major neutral and ionic species according to a standard 1D SIC run.The aim of the reduction was to find the smallest sub-mechanism which has a maximum root-mean-square error ( RMS ) of important species at each of the four altitudes less than 5% (i.e. RMS <0.05).Note that the 5% accuracy is guaranteed only for the selected 1D conditions (altitudes = 60, 70, 80, 90 km) but not for any other conditions.The 1D model was run for 10 days and local errors were calculated at 100 points distributed logarithmically between 1 s and 10 days (864,000 s), with an increase of approximately 15% between adjacent times.In order to set the proton flux to peak, the reduction was done at the peak time of the storm (0615 UT, 29 th October 2003).

3D modelling using WACCM
WACCM is a comprehensive numerical model extending vertically from the ground up to the lower thermosphere (~140 km) and is part of the NCAR Community Earth System Model (CESM) (Hurrell et al., 2013).Here we use the specified dynamics (SD) version of WACCM 4 (Marsh et al., 2013), which has 88 pressure levels from the surface to 5.96×10 -6 hPa and a horizontal resolution of 1.9 o × 2.5 o (latitude × longitude).The model contains all the important processes in the atmosphere: chemistry, radiative transfer, auroral processes, non-local thermodynamic equilibrium, ion drag, SPE ionization rate from 1963 to 2016 and the molecular diffusion of the constituents.In the SD version the model is forced with ECMWF meteorology re-analysis (Dee et al., 2011) from surface to 50 km.1% of the meteorological conditions (temperature, winds, surface pressure, specific humidity, surface wind stress, latent, sensible heat flux) were combined with the WACCM fields below 50 km at every model dynamics time step.This nudging factor then reduces linearly from 1% to 0% between 50 km and 60 km.Above 60 km there is no nudging to the re-analysis fields and the model is free-running.The EPP caused NO x production was determined from the ionization rate (I) that is proportional to the energy deposition rate, and the NO x production rate is expressed as 1.25×I.This is the standard parameterization in WACCM using Jackman's method: According to (Jackman et al., 2005) an ion pair produces 1.25 N atoms with branching ratios of 0.55 N( 4 S) and 0.70 N( 2 D).
Then the nitrogen atoms will go through a complex ion chemistry that leads to the formation of NOx species..For The HO x there was no similar parameterization involved species are produced via a complicated ion chemistry scheme that is based on (Solomon et al., 1981).
Both the full SIC model and its reduced version, rSIC, were included in WACCM to produce the WACCM-SIC and WACCM-rSIC models, respectively.WACCM-rSIC included the WACCM neutral reactions and the rSIC module which has a similar number of reactions ( 181 Cl -(HCl) and NO 3 -(HNO 3 )).In the case of WACCM-rSIC the three major positive and negative ions (NO + , H + (H 2 O) 3 , -and NO 3 -) were selected (based on the results from the 1D model run for January 2005), and only the two-body ion-ion reactions of these species with the non-redundant ion species (of opposite charge) were included.For WACCM-SIC, the two-body ion-ion reactions of these three major positive and three major negative ions with all ions of opposite charge were included.This is a good approximation as this still allows the ion-ion channels to contribute to the net ion balance but ignores those channels where the contributing ion concentrations do not dominate.Three-body ion-ion recombination reactions are many orders of magnitude slower than the two-body processes in the D region, and hence were not included in both the full and the reduced 3D models.
After a 5-year spin up, all four models -WACCM, WACCM-SIC, WACCM-rSIC and WACCM-Dwere run for 30 days from 1 st January 2005, a period encompassing a significant SPE which reached a maximum on 17 th January 2005 at 17:50 UT.Long spin up period was used in order to definitely allow enough time to reach steady-state and also to count for long lived species such as N 2 O, CH 4 .
The January 2005 event (15 th January 2005 -20 th January 2005) was chosen in order to test the accuracy of the reduced model for a different scenario to the Halloween storm.The reason for choosing this event was to be consistent with our earlier study using WACCM-D (Verronen et al., 2015).

Results and discussion
The concentration tolerance for the important species was set to 5%, meaning that the concentrations of the important species in the reduced model are reproduced to within 5% of the full model.The change of root-mean-square (rms) error as a function of number of species is shown in Figure 1, where it is clear that the error drops initially rapidly with increasing number of species and then continues to decrease more slowly.The species were added in order of decreasing error.Until the reaction pathways are built up completely in the reduced mechanism, various feedback effects are not necessarily manifested in the concentration profiles of the important species.The intermediate reduced mechanism will therefore have one or more inactive truncated pathways.A sudden decrease in the rms error is observed after a long stagnation when those species are added which close the pathways; their early inclusion would not give a smaller reduced mechanism, as then other species would be missing from the pathway.The jump in the numbers of necessary species from 60 km to 70 km is connected to the increased concentrations of NO + and O 2 + -these ions are involved in several ion-molecule reactions.The sets of important and necessary species are tabulated in Table 1.These were selected on the basis of which species are expected to be chemically directly associated with NO, NO 2 , HNO 3 and HO 2 ; the five most abundant positive and negative ions were essential to be included.Ions were also selected in the initial set of important species.When a neutral-only list of important species was tried the reduction never reached the required accuracy which indicated that ions are crucial part for the balance of the major neutrals via neutral-ion interconnection.Bold-face-font species indicate the initial set of important species, while the others are the necessary species selected through the reduction process.As electrons are not important charge carriers in the D region, they were not considered as important species.The reduced 1D model decreases the computational time by 25% compared to the full model, while the number of reactions dropped from 306 307 to 181 and the number of species from 129 to 81.However, more importantly, the results of the reduction process makes it possible to identify the reactions of the important and necessary species that affect the atmospheric NO x and HO x concentrations, and can therefore impact on stratospheric ozone levels.

1D modelling
Figure 2 compares the major neutral and ion profiles produced by the SIC and rSIC models for the SPE conditions of January 2005.As expected, the concentrations of important neutral and ionic species are reproduced very well by the reduced model, within the height range (60 -90 km) where the reduction was carried out.Figure 2 also shows the relative errors of the mechanism reduction on the concentration profiles.Note that the differences can only be seen in the right hand panels of Figure 2. It is clear that the errors remain within the 5% tolerance, except for HO 2 below 30 km and for O 2 -above 115 km, but these are outside the D region.Furthermore, the neutral profiles are reproduced within a factor of 52 under any conditions, even above 80 km.This is not surprising as the chemical scheme in the lower thermosphere is much simpler, being governed by NO + , O 2 + , O + and electrons (Danilov, 1994).
All photodissociation and photoionization reactions are found to be important, and are therefore included in the 181 reactions in the reduced model.27% of the positive ion recombination reactions could be eliminated, which are those involving recombination of NO + (N 2 ), NO + (CO 2 ) and their mono-and di-hydrated complex clusters.The 7 and 8 H 2 O-containing proton hydrates are also redundant species.All the other reactions involving electrons, including photodetachment and reactive electron detachment are important.This will be discussed later in more detail.The reduction was clearly most effective for the positive ions with 55 out of the 124 positive ion reactions (i.e.44%) being redundant, including: the reactions of CH 3 CN and all its clusters; all the proton hydrate + N 2 O 5 reactions; and the reactions of proton-hydrate-HNO 3 complex clusters with N 2 O 5 and/or H 2 O.The majority of the redundant species are the hydrated clusters, and only a few N 2 and CO 2 clusters are included, such as the reactions of NO + (N 2 ) and NO + (CO 2 ).
For the negative ions, 54 out of the total of 124 reactions (i.e.43%) were redundant, which is similar to the subset of redundant positive ion reactions.Similar to the positive ion results, the redundant list mostly arises from the reactions of hydrate clusters.From the reduced set it can be concluded that in general the majority of the positive ion reactions are the reactions of NO + or O 2 + and their hydrated clusters.The reactions of proton-hydrate clusters up to n = 6 also turn out to be necessary under all conditions, apart from the H + (H 2 O) n + N 2 O 5 reactions; the hepta-and octa-hydrates are both redundant species.The mechanism reduction in the vicinity of the boundary between the mesosphere and lower thermosphere is more efficient, with 60% of the reactions found to be redundant.This is not surprising as lower thermospheric chemistry is governed chiefly by NO + and O 2 + , and the scheme simplifies as the negative ion concentrations (principally O 2 -, CO 3 -and NO 3 -) decrease by many orders of magnitude from 80 km to 90 km.The reactions of HCO 3 -with neutrals are redundant under all conditions.OH -generation and consumption are important only at ~60 km, and above this height OH -becomes redundant.The reactions of CO 4 -are redundant already between 70 km and 80 km, while the NO 2 -reactions are redundant between 80 km and 90 km.The total number of reactions in the reduced model is 181.Testing the reduced mechanism with the 3D model is described below.

3D WACCM modelling
The full and reduced SIC models were coupled into WACCM to produce WACCM-SIC and WACCM-rSIC, respectively.
The standard WACCM model contains only the five major positive ions (N + , N 2 + , O + , O 2 + and NO + ) and has no negative ions.WACCM-D is described in our recent paper (Verronen et al., 2016).The WACCM-D reactions are shown in bold typeface in the SM.From this it is clear that the majority of reactions in WACCM-rSIC and WACCM-D are the same.
However, WACCM-D omits the reactions of the larger proton hydrates and their clusters ).These four WACCM-based models (WACCM, WACCM-D, WACCM-SIC and WACCM-rSIC) were tested for the 17 th January 2005 SPE, a moderate event with peak flux of 5040 (pfu>10MeV) at 17:50 UTC.For auroral electrons and SPEs, the vertical profile of the energy deposition is inferred from electron and proton fluxes, observed in low earth orbit and in geostationary orbit, respectively (Holt, 2013).
We investigated the effect of the SPE on O 3 and the major NO x and HO x species.The results are shown in Figure 3-7, for northern hemisphere polar latitudes (60 -90N) during January 2005.Figure 3a shows the calculated polar O 3 mixing ratio from the four WACCM based models for January 2005; the percentage differences from WACCM-SIC are shown in Figure 3b.It is important to note that the 5% accuracy was guaranteed only for the selected 1D conditions.Therefore the deviation can be quite different at other conditions which means that the > 5% accuracy with the 3D model is not surprising.On the other hand, the different treatments of the ion-ion reactions in WACCM-SIC and WACCM-rSIC could also increase the deviationsie.only the reactions of the major positive and negative ions with the non-redundant species were selected for the WACCM-rSIC, while for WACCM-SIC all the major positive and negative ion reactions were considered.Inspection of this figure shows generally very good agreement (average difference <10%, maximum difference 30%) between WACCM-rSIC and WACCM-SIC.There are somewhat larger differences between WACCM-D and WACCM-SIC, with an overestimate of the O 3 concentration by up to 100% between 65 and 70 km.This is due to the lack of proton hydrate and O 2 + cluster reactions in WACCM-D.Note that the standard WACCM fails to simulate the O 3 depletion between 55 and 70 km during the SPE (January 16 th -23 rd ), which illustrates the effect of not including positive cluster ions and all negative ions in WACCM.(Jackman et al., 2011) reported the concentrations of OH, HO 2 , HNO 3 and O 3 measured by the Microwave Limb Sounder (MLS) satellite instrument, and NO x data measured by the Advanced Composition Explorer (ACE).In the lower mesosphere (55 -70 km), O 3 decreased by 15-60% (with the largest change at ~12:00 UT on 17 th January).Figure 3a shows that this behaviour is captured well by WACCM-rSIC, WACCM-D and WACCM-SIC, which all exhibit a decrease of 20-60% in only about 12 hours.
The polar NO, NO 2 , OH, and HO 2 vertical profiles are shown in Figures 4-7.There is a marked increase in NO by a factor of 5 -10 between 45 and 70 km, starting directly after the beginning of the SPE (02:10 UT on 16 th January) and lasting for more than 15 days (Figure 4a).The NO then recovers back very slowly (~30-day timescale) to the original concentration.
The sudden increase of NO concentration after the SPE is picked up even by standard WACCM, indicating that this NO production is largely governed by the chemistry of the five major positive ions; however, standard WACCM predicts the NO enhancement down to 63 km, in contrast to the three SIC models where NO enhancement only occurs down to 70 km.This is most likely due to the significant role of negative ions in the 60-70 km region: the major negative ions, O -, CO 3 -and ClO -, convert NO to NO 2 or NO 2 -(reactions NIR40, NIR57, NIR122 and NIR123 in the SM).
Figure 5a shows that the NO 2 density increases by a factor of ~3 after the beginning of the SPE.The effect is also seen in the standard WACCM run, but the effect is overestimated by a factor of 2.5 compared to WACCM-SIC.This is due to the lack of anion chemistry: NO 2 is converted to NO 2 -via reaction with O -, O 2 -and Cl -(NIR4, NIR17, NIR113), or to NO 3 -via reaction with CO 3 -(NIR58 in the SM). Figure 5b shows that WACCM-D overestimates NO 2 on average by 40% compared to WACCM-SIC and WACCM-rSIC.This appears to be due to the absence of N 2 clusters of proton hydrates and of the O 2 + ion in WACCM-D.The relatively small difference between WACCM-rSIC and WACCM-SIC arises from the absence of the reactions of hydrated O 2 -ions (NIR(28) and NIR(31) in the SM).In general, the pre-SPE NO and NO 2 concentrations from WACCM-D are more similar to WACCM than to WACCM-rSIC.Also, the post-SPE recovery of NO is more similar to the WACCM results, while the NO 2 recovery shows a mixed behaviour compared to WACCM-rSIC and WACCM: the increase rate at the SPE is around half way compared to the WACCM-rSIC and WACCM results.These is explained by the different NO x chemistry in WACCM-D and WACCM-rSIC: as it can be seen in the SM, a the important recombination reactions of NO + clusters forming NO, namely RPE(8) and RPE(9) are missing in WACCM-rSIC.Although these do not include direct NO 2 formation, the increase in NO concentration will affect the NO 2 concentrations in both WACCM-D and WACCM-rSIC via the negative ion reactions of NIR(40), NIR(57).
Figure 6a shows that WACCM-SIC, -rSIC and -D all predict OH increases by a factor of ~5 between 53 and 70 km during the SPE, and a rapid recovery to quiet time levels after the SPE ends on January 23 rd .In contrast, standard WACCM exhibits a modest enhancement of OH, and only above 63 km. Figure 7a shows HO 2 enhancements up to a factor of ~3 during the SPE, extending all the way down to 40 km -apart from standard WACCM.Again, the lack of anion chemistry seems to be responsible, since these HO x species are produced both by electron detachment (EDA(9), EDA(12) and EDA( 13 and 80 km during the maximum of the SPE, and an even larger increase of up to 150% between 60 and 70 km (Jackman et al., 2011).In comparison, WACCM-SIC and -rSIC predict increases of 75% and 150% in these respective height ranges, which are therefore in sensible agreement with the satellite observations.(Jackman et al., 2011) concluded that this large increase in HO 2 was mostly responsible for the catalytic destruction of O 3 in the mesosphere.
Figure 8a shows the HNO 3 density predicted by the four models.Standard WACCM is completely unable to reproduce the HNO 3 concentration enhancement predicted by WACCM-SIC, which underlines the importance of negative ion chemistry.
The differences between WACCM-rSIC and WACCM-D are generally smaller than 50% and the effect of the SPE on HNO 3 is similar in both reduced models.Compared to WACCM-SIC, WACCM-rSIC overestimates HNO 3 by over a factor of 2 between 55 and 62 km during the SPE, while WACCM-D also overestimates HNO 3 but over a wide altitude range for a Therefore, the extra loss of HNO 3 via (R1) and (R2) are overcompensated by these HNO 3 increasing channels.(Verronen et al., 2011) confirmed earlier work (Aikin, 1997;Verronen et al., 2008) which showed that the most important HNO 3 -forming reactions in the mesosphere are the H + (H 2 O) n + NO 3 -(HNO 3 ) m → (m+1)HNO 3 + n H 2 O ion-ion recombination reactions.We found all these reactions to be necessary.However, the positive ion channels of direct HNO 3 formation are all redundant (i.e.(PIR97)-( PIR101) and ( PIR102)-(PIR106) in the SM): The NO 3 -(H 2 O) n + N 2 O 5 → HNO 3 channel is also inefficient under all conditions.The NO 3 -ion or its hydrates are responsible for HNO 3 formation only via their reactions with proton hydrates, while CO Therefore, in WACCM-rSIC only the negative ion reactions affect the HNO 3 balance, as no positive ions contribute to its production or its loss.In contrast, the concentrations of O 3 , NO x and HO x are influenced by both positive and negative ion reactions which remain in the reduced model (see the SM for further details).It should be noted that the identification in the present study of unnecessary reactions of N 2 O 5 allows a sub-mechanism to be extracted that still reproduces the mesospheric O 3 , HNO 3 , HO x and NO x profiles with reasonable accuracy.
A detailed analysis of the ionic reaction sequences and resulting changes in neutral species has been provided by (Verronen and Lehmann, 2013), who concluded that positive ion chemistry mainly leads to OH and NO formation, while negative ion chemistry causes HNO 3 enhancement via NO, NO 2 , N 2 O 5 → HNO 3 conversion.This is partly in accordance with our finding that the HNO 3 concentration is determined by negative ion reactions.However, we found that the OH concentration is also influenced by negative ion reactions of O -, OH -and CO 3 -(reactions PDE(4), NIR(3), NIR(7), NIR(11), NIR(47), NIR(48), NIR(52), NIR(53), NIR(59), NIR(64) and NIR(73) in the SM).According to Verronen and Lehmann, mesospheric NO production is dominated by positive ion reactions; here we conclude that NO is consumed mostly by negative ion reactions.

Conclusions
In WACCM-rSIC runs 25% faster than WACCM-SIC, which is a reasonable improvement.Moreover, the simulation time with WACCM-rSIC is only 90% more than with the standard version of WACCM, which contains no negative ions or positive cluster ions and is clearly inadequate for simulating the impacts of energetic particle precipitation on the chemistry of the middle atmosphere.Note that the rSIC model contains temperature (and pressure) dependent rate coefficients (where available), therefore the chemical model can be used over a wide range of conditions.
this paper we have described the development of a whole D region chemistry based on the Sodankylä Ion and Neutral Chemistry (SIC) coupled into WACCM.Systematic mechanism reduction using the Simulation Error Minimization Connectivity method produces a D region ion chemistry described by 181 ion-molecule reactions, a 41% reduction from the full SIC chemistry.Our earlier WACCM-D version has 196 reactions; the majority of the reactions in WACCM-rSIC and WACCM-D are the same.However, WACCM-D does not include reactions of the larger proton hydrates and their clusters, as well as the reactions of O 2 + clusters, and inclusion of these reactions gives better agreement with the O 3 predicted by the full WACCM-SIC chemistry.Before and after the solar proton event, NO and NO 2 from WACCM-rSIC agree much better with results from WACCM-SIC than the results of WACCM-D, which is due to the absence of N 2 clusters of proton hydrates and of the O 2 + ion in WACCM-D.

Figure 1 .Figure 2 .
Figure 1.Root-mean-square (rms) error of all species as a function of number of species in the reduced mechanism for the four selected altitudes (60, 70, 80 and 90 km). 5

Table 2 .
Photodissociation and photoionization reactions of the SIC model