D-region ion – neutral coupled chemistry ( Sodankylä Ion Chemistry , SIC ) within the Whole Atmosphere Community Climate Model ( WACCM 4 ) – WACCM-SIC and WACCM-rSIC

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 Chemistry (SIC) model, a one-dimensional model containing 307 ion–neutral and ion recombination, 16 photodissociation and 7 photoionization reactions of neutral 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 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 (HNO3, O3, H2O2, NO, NO2, HO2, OH, N2O5) and ions (O+2 , O + 4 , NO , NO(H2O), O+2 (H2O), H (H2O), H(H2O)2, H(H2O)3, H(H2O)4, O−3 , NO − 2 , O , O2, OH, O−2 (H2O), O − 2 (H2O)2, O − 4 , CO−3 , CO − 3 (H2O), CO − 4 , HCO − 3 , NO − 2 , NO − 3 , NO − 3 (H2O), NO−3 (H2O)2, NO − 3 (HNO3), NO − 3 (HNO3)2, Cl , ClO), which agree with the full SIC mechanism within a 5 % tolerance. Four 3-D 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 WACCMrSIC (standard WACCM with a reduction of SIC chemistry using the SEM-CM method). The 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 ion–molecule reactions that affect the partitioning of odd nitrogen (NOx), odd hydrogen (HOx) and O3 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; in the case of electrons, from the radiation belts around the Published by Copernicus Publications on behalf of the European Geosciences Union.
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.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., 2007Marsh 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 leading kinetic models of D-region chemistry is the Sodankylä Ion 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 307 ion-molecule and also 2254 ionion recombination reactions.The model is a 1-D model with simple vertical transport (molecular and eddy diffusion) to balance the computational cost.Putting such a large number of additional reactions into a 3-D 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 schemei.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 the 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 mediumintensity 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 307 ionneutral reactions in total of all the important and necessary species.These are listed in the Supplement, 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 Supplement.

Description of the 1-D 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 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).The SIC model inputs ionization rates as a function of time and pressure, choosing 3 h time resolution.The Mass Spectrometer and Incoherent Scatter Radar Extended (MSISE-90) model (Picone et al., 2002) is used for climatological average altitude profiles of N 2 , O 2 and temperature for the background atmosphere, and midlatitude concentrations of N 2 O, H 2 , HNO 2 , HCl, Cl, ClO, CH 3 , CH 4 , CH 2 O and CO at altitudes of 10,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 species -that is, the non-redundant reactions.The set of redundant species takes 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: dc dt = f (c, k(T , p), J (F (λ), T , p)) . (1) 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 (307).Note that the photoionization and photodissociation reactions are not listed in the Supplement 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 ( Jij = ∂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 kth 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 an initially present species or whether 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"); 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 are 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 were found to be very efficient.These errors were defined as δ RMS = rms Here δ ij k is a local dimensionless error measuring the deviation between the concentrations of important species j at time point k in scenario i modelled by the full (c full ij (t k )) and the reduced mechanisms (c red ij (t k )).A mixed type of local error defined by Nagy and Turanyi (2009) was used: This mixed error behaves as a relative error close to the maximum concentration of species j in scenario i (c full ij ,MAX ) 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 06:15 UT on 29 October 2003 with a proton flux of 29 500 (pfu > 10 MeV) (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 1-D 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 1-D conditions (altitudes = 60, 70, 80, 90 km), not for any other conditions.The 1-D 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 (06:15 UT, 29 October 2003).

3-D 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 • × 2.5 • (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 European Centre for Medium-Range Weather Forecasts (ECMWF) meteorology re-analysis (Dee et al., 2011) from the surface to 50 km.One per cent of the meteorological conditions (temperature, winds, surface pressure, specific humidity, surface wind stress, latent, sensible heat flux) was 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 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 go through a complex ion chemistry that leads to the formation of NO x species.The HO x 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) to , but they differ with the reactions of O + 2 clusters and larger proton hydrates (see the Supplement).As ion-ion reactions were not included in the reduction process, they were considered in a different way for WACCM-D, WACCM-SIC and WACCM-rSIC.WACCM-D has 112 two-body and 14 three-body reactions.The two-body reactions involves the In the case of WACCM-rSIC the three major positive and negative ions (NO 2 , CO − 3 and NO − 3 ) were selected (based on the results from the 1-D 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 3-D models.

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 Fig. 1, where it is clear that the error initially drops 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 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 a 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 1-D model decreases the computational time by 25 % compared to the full model, while the number of reactions dropped from 307 to 181 and the number of species from 129 to 81.However, more importantly, the results of the reduction process make 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.

1-D 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 righthand panels of Fig. 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 5 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.Twenty-seven per cent 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 seven and eight 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 rewww.geosci-model-dev.net/9/3123/2016/Geosci.Model Dev., 9, 3123-3136, 2016 Geosci.Model Dev., 9, 3123-3136, 2016 www.geosci-model-dev.net/9/3123/2016/ Table 2. Photodissociation and photoionization reactions of the SIC model.
actions; 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 clus-ters.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 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 and 80 km, while the NO − 2 reactions are redundant between 80 and 90 km.The total number of reactions in the reduced model is 181.Testing the reduced mechanism with the 3-D model is described below.

3-D 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 Supplement.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 clus-  a moderate event with peak flux of 5040 (pfu > 10 MeV) 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 et al., 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 Figs.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 Fig. 3b.It is important to note that the 5 % accuracy was guaranteed only for the selected 1-D conditions.Therefore the deviation can be quite different at other conditions, which means that the > 5 % accuracy with the 3-D 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 deviations -i.e.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 (16-23 January), 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 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 h.
The polar NO, NO 2 , OH and HO 2 vertical profiles are shown in Figs.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 January) and lasting for more than 15 days (Fig. 4a).The NO then recovers very slowly (∼ 30-day timescale) to the original concentration.The sudden increase of NO concentration after the SPE is picked up even by the standard WACCM, indicating that this NO production is largely governed by the chemistry of the five major positive ions; however, the 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 − 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 2 ions (NIR(28) and NIR(31) in the Supplement).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 halfway compared to the WACCM-rSIC and WACCM results.This is explained by the different NO x chemistry in WACCM-D and WACCM-rSIC: as can be seen in the Supplement, the 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) and NIR(57).
Figure 6a shows that WACCM-SIC, WACCM-rSIC and WACCM-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 23 January.In contrast, the 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 the 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 anion-molecule reactions (NIR(3), NIR(6), NIR(7), NIR(11), NIR(20), NIR(23), NIR(34), NIR(47), NIR(48), NIR(50), NIR(52), NIR(53), NIR(56), NIR(59), NIR(64), NIR(73), NIR(77) and NIR(80)).Figures 6b and 7b show that WACCM-rSIC generally produces better agreement with WACCM-SIC than WACCM-D.MLS observed a doubling of HO 2 between 70 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 WACCM-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.The 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 greater duration.The only HNO 3consuming reactions that are eliminated in WACCM-rSIC are two fast HNO 3 -cluster-forming Reactions (R1) (NIR99) WACCM-SIC and WACCM-rSIC differ only in the ion chemistry.This indicates that the difference in the concentrations can arise only from this.Although Reactions (R1) and (R2) are present in WACCM-D, and therefore one would expect a smaller response of HNO 3 during the SPE, from the Supplement it is clear that there are a couple of indirect chan- 3 (H 2 O) raises the HNO 3 concentration via NIR(97), which reaction also increases the NO − 3 (HNO 3 ) concentration that raises the HNO 3 concentration via NIR(101).Therefore, the extra loss of HNO 3 via Reactions (R1) and (R2) is overcompensated for 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 ionion 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 Supplement): (R4) 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 − 3 , O − , O − 2 and Cl − are responsible for HNO 3 removal via Reactions (R5)-(R8) ((NIR64), (NIR11), (NIR23), (NIR115)).
The NO − 2 → NO − 3 conversion also removes HNO 3 by converting it into HNO 2 via Reaction (R9) (NIR86) below 80 km: 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 Supplement 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 and 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 Supplement).According to Verronen and Lehmann, mesospheric NO production is domi-nated by positive ion reactions; here we conclude that NO is consumed mostly by negative ion reactions.

Conclusions
In this paper we have described the development of a whole D-region chemistry based on 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, nor 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.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.
The combined WACCM-SIC and the reduced version, WACCM-rSIC, as well as all model simulation results are available upon request to J. M. C. Plane, while the WACCM-D and the SIC model are available upon request to P. T. Verronen.
The Supplement related to this article is available online at doi:10.5194/gmd-9-3123-2016-supplement.

Figure 1 .
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).

Figure 2 .
Figure 2. Atmospheric concentration profiles and relative differences calculated using the full SIC model and reduced (rSIC) model for 17 January 2005 at 17:50 UT.In the left-hand panels the solid lines show the concentrations calculated by the full SIC model, while the symbols refer to the reduced rSIC model.The right-hand panels show the percentage difference between the rSIC and SIC models: (a) concentrations and (b) percentage differences of HNO 3 , NO, NO 2 , O 3 and H 2 O 2 ; (c) concentrations and (d) percentage differences of OH and HO 2 ; (e) concentrations and (f) percentage differences of NO + and O + 2 ; and (g) concentrations and (h) percentage differences of CO − 3 , NO − 3 and O − 2 .