Statistical response of middle atmosphere composition to solar proton events in WACCM-D simulations: the importance of lower ionospheric chemistry

Atmospheric effects of solar proton events (SPEs) have been studied for decades, because their drastic impact can be used to test our understanding of upper stratospheric and mesospheric chemistry in the polar cap regions. For example, odd hydrogen and odd nitrogen are produced during SPEs, which leads to depletion of ozone in catalytic reactions, such that the effects are easily observed from satellites during the strongest events. Until recently, the complexity of the ion chemistry in the lower ionosphere (i.e., in the D region) has restricted global models to simplified parameterizations of chemical impacts induced by energetic particle precipitation (EPP). Because of this restriction, global models have been unable to correctly reproduce some important effects, such as the increase in mesospheric HNO3 or the changes in chlorine species. Here we use simulations from the WACCM-D model, a variant of the Whole Atmosphere Community Climate Model, to study the statistical response of the atmosphere to the 66 strongest SPEs which occurred in the years 1989–2012. Our model includes a set of D-region ion chemistry, designed for a detailed representation of the atmospheric effects of SPEs and EPP in general. We use superposed epoch analysis to study changes in O3, HOx (OH + HO2), Clx (Cl + ClO), HNO3, NOx (NO + NO2) and H2O. Compared to the standard WACCM which uses an ion chemistry parameterization, WACCM-D produces a larger response in O3 and NOx and a weaker response in HOx and introduces changes in HNO3 and Clx . These differences between WACCM and WACCM-D highlight the importance of including ion chemistry reactions in models used to study EPP.

Abstract. Atmospheric effects of solar proton events (SPEs) have been studied for decades, because their drastic impact can be used to test our understanding of upper stratospheric and mesospheric chemistry in the polar cap regions. For example, odd hydrogen and odd nitrogen are produced during SPEs, which leads to depletion of ozone in catalytic reactions, such that the effects are easily observed from satellites during the strongest events. Until recently, the complexity of the ion chemistry in the lower ionosphere (i.e., in the D region) has restricted global models to simplified parameterizations of chemical impacts induced by energetic particle precipitation (EPP). Because of this restriction, global models have been unable to correctly reproduce some important effects, such as the increase in mesospheric HNO 3 or the changes in chlorine species. Here we use simulations from the WACCM-D model, a variant of the Whole Atmosphere Community Climate Model, to study the statistical response of the atmosphere to the 66 strongest SPEs which occurred in the years 1989-2012. Our model includes a set of D-region ion chemistry, designed for a detailed representation of the atmospheric effects of SPEs and EPP in general. We use superposed epoch analysis to study changes in O 3 , HO x (OH + HO 2 ), Cl x (Cl + ClO), HNO 3 , NO x (NO + NO 2 ) and H 2 O. Compared to the standard WACCM which uses an ion chemistry parameterization, WACCM-D produces a larger response in O 3 and NO x and a weaker response in HO x and introduces changes in HNO 3 and Cl x . These differences between WACCM and WACCM-D highlight the importance of including ion chemistry reactions in models used to study EPP.

Introduction
Solar proton events (SPEs) are observed on Earth when highenergy protons accelerated in the sun's magnetic field during a solar coronal mass ejection strike Earth. These protons generally have energies in the range of tens or hundreds of megaelectronvolts per nucleon and due to their high energies are able to penetrate deep into the atmosphere at geomagnetic latitudes polewards of ≈ 60 • . This causes ionization and dissociation (mainly of the most abundant species N 2 and O 2 ) in the altitude range of roughly 30-90 km. Solar proton events typically last from a few hours to few days and can occur anytime during the 11-year solar cycle, although they are more common during solar maximum.
Ionization, excitation and dissociation caused by the SPEs, and energetic particle precipitation (EPP) in general, have a significant influence on neutral composition through ionneutral chemistry (e.g., Verronen and Lehmann, 2013). Sev-N. Kalakoski: SPE impact on atmospheric composition eral studies have investigated the depletion of ozone in the polar mesosphere and upper stratosphere resulting from the increased production of odd hydrogen and odd nitrogen (e.g., Jackman et al., 2001;Seppälä et al., 2004;López-Puertas et al., 2005;. Ozone depletion can result in changes in temperatures during the largest SPEs , and there is evidence that SPEs, as part of EPP, can modulate polar winter dynamics on decadal timescales Semeniuk et al., 2011;Rozanov et al., 2012;Calisto et al., 2013). Recent results have also shown that SPEs play a role in the upper stratospheric variability of ozone and have to be considered when studying the expected recovery of the ozone layer (Stone et al., 2018). The state of the art of the impact of energetic particle precipitation, in particular SPE and geomagnetic forcing, has also been summarized in the recent WMO assessment (WMO, 2018, chap. 4.3.5.1, 4-25, and references therein).
Simulation results have suggested that the decrease in the polar total ozone column would be of the order of 1 %-2 % a few months after large SPEs (Jackman et al., 2014). However, local depletion has been reported to reach ≈ 10 % below the ozone layer peak at 50-100 hPa, based on a statistical analysis of almost 200 SPEs using ozone soundings (Denton et al., 2018). As the contribution of >300 MeV protons to direct ozone loss in the lower stratosphere would likely be negligible due to the relatively small fluxes at such high energies , any observed ozone loss would likely be the result of indirect effects. Ozone may also increase in the lower stratosphere due to the enhanced NO x interfering with chlorine-driven catalytic ozone loss (Jackman et al., 2008).
The lower ionospheric, the D-region, chemistry, which connects SPE ionization to changes in neutral species, is rather complex compared to the five-ion chemistry that is adequate in the thermosphere. A special feature of the D-region is the presence of negative ions which are formed when electrons attach to neutral species. Another feature is the large number of different cluster ions, both positive and negative, some of which are among the most abundant ions in the region. Due to the large number of D-region ions and ionic reactions, atmospheric models have typically included the SPE, or EPP, effects using simple parameterizations of HO x and NO x production. Funke et al. (2011) discussed shortcomings related to these parameterization schemes in a comparison study between observations and simulations of the October-November 2003 SPE. Among the outstanding issues in simulations have been the lack of nitric acid HNO 3 increase (see also Jackman et al., 2008) and the lack of chlorine activation. The HNO 3 production is driven through reactive nitrogen redistribution by negative ion chemistry (Verronen and Lehmann, 2013), and simulations could be improved either by improved parameterization , or by considering the relevant ion chemistry explicitly (Verronen et al., 2008;Verronen et al., 2011;Andersson et al., 2016). The differences between modeled and observed response of chlorine species to SPE reported by Jackman et al. (2008) can be explained by the conversion of HCl to active chlorine species (Cl, ClO, HOCl) by the ion chemical reactions (Winkler et al., 2009(Winkler et al., , 2011. In the lower ionosphere, the ion chemistry reactions are expected to lead to the depletion of water vapor during SPEs. For moderate-sized SPEs, this effect is small compared to the icy particle sublimation governed by the changes in temperature (von Savigny et al., 2007;Winkler et al., 2012).
The Whole Atmosphere Community Climate Model (WACCM) is a global chemistry-climate model and forms the atmospheric part of the Community Earth System Model (CESM). Recently, a variant called WACCM-D was developed for detailed simulations of D-region ion chemistry and EPP atmospheric effects . In a case study of the January 2005 SPE, the consideration of ion chemistry in WACCM-D led to an improved SPE response in HNO 3 , Cl x , NO x and HO x . Since then, WACCM-D has been used to study mesospheric nitric acid and cluster ion composition during electron precipitation events (Orsolini et al., 2018) and magnetic latitude dependency of SPE ionospheric impact (Heino et al., 2019).
Here, we take a statistical approach to look at the response of middle atmosphere chemistry to a number of SPEs of various intensities. As the SPE effects are largely known from previous work, we focus on the improvement provided by additional ion chemistry reactions included in the WACCM-D chemistry. In addition to the traditionally analyzed species, such as O 3 , HO x , NO x and HNO 3 , we also investigate water vapor and active chlorine species which have been less studied, especially with global models. While several of the largest SPEs included in the analysis have previously been studied individually, the statistical approach used here allows for inclusion of a number of moderate-sized events under various background atmosphere and illumination conditions. This approach also allows for the identification of climatological effects above natural variability. As the analysis includes SPEs of different sizes occurring during different seasons, a statistical approach is most useful for the study of temporal and spatial extent, rather than magnitude of the response.
2 Modeling and analysis methods

WACCM-D simulations
WACCM is a global circulation model with fully coupled chemistry and dynamics extending from the surface to 6 × 10 −6 hPa (≈ 140 km). The horizontal resolution used is 1.9 • latitude by 2.5 • longitude. A description of the model physics in the MLT (mesosphere-lower thermosphere) as well as the simulations of dynamical and chemical response to radiative and geomagnetic forcing during solar maximum and min-imum are described by . Marsh et al. (2013) present an overview of the model climate and describes the climate and the variability in a long-term simulation using version 4 of the WACCM. The standard chemistry package of WACCM includes photochemistry of 59 neutral species for the whole altitude range, and five ion species O + , NO + , O + 2 , N + 2 and N + for the lower thermosphere, as well as excited species N( 2 D), O( 1 D), O 2 ( 1 ) and O 2 ( 1 ). For SPE effects, and energetic particle precipitation in general, HO x and NO x production is parameterized. For the SPE NO x effect, it is assumed that 1.25 N atoms are produced per ion pair with branching ratio of 0.55/0.7 for N( 4 S)/N( 2 D) Porter et al., 1976). This parameterization depends on a fundamental assumption of a fixed N 2 /O 2 ratio, and it has been shown to underestimate NO x production above about 65 km (Nieder et al., 2014;Andersson et al., 2016). For the SPE HO x effect, based on the work of Solomon et al. (1981), a maximum of two molecules are produced per ion pair in the lower mesosphere and below while in the upper mesosphere the production gradually decreases to zero with increasing altitude.
WACCM-D is a variant of WACCM in which the standard parameterizations of HO x and NO x production are replaced by a set of lower ionospheric photochemistry, with the aim of reproducing better the observed effects of EPP on the mesosphere and upper stratosphere neutral composition. The ion chemistry set was selected based on analysis of the latest knowledge of chemical reactions in the ionospheric Dregion and their effects on the neutral atmosphere (Verronen and Lehmann, 2013). The set consists of 307 reactions of 20 positive ions and 21 negative ions, including cluster ions such as H + (H 2 O) and NO − 3 (HNO 3 ), which are important for HO x and HNO 3 production, for example. The details on WACCM-D ion chemistry as well as its lower ionospheric evaluation were presented by Verronen et al. (2016), while the improvement in the atmospheric response during the January 2005 SPE was demonstrated in comparisons with satellite observations . WACCM-D is now included as an official predefined component set (compset) since the CESM 2.0 release in June 2018.
In this study, we use WACCM version 4 simulations with the specified dynamics configuration (SD-WACCM) which is forced with meteorological fields (temperature, horizontal winds and surface pressure) from NASA's Modern-Era Retrospective Analysis for Research and Applications (MERRA) (Rienecker et al., 2011;Lamarque et al., 2012) at every dynamics time step below about 50 km; above this, the model is free running (88 levels in total). It should be noted that even above 50 km model dynamics are strongly modulated by winds and wave fluxes from lower levels. This means that internal variability of dynamics in SD-WACCM is small.
For energetic particle precipitation, the simulations include forcing from auroral electrons (E < 10 keV), solar protons and galactic cosmic rays. Medium-energy and high-energy electrons (E > 10 keV) were excluded from the simulations. Two model simulations were made: (1) a reference run using a standard SD-WACCM compset (referred to as REF) and (2) a run with D-region ion chemistry (referred to as WD). Both simulations covered the years 1989 to 2012. We use daily mean SPE ionization rates based on GOES (Geostationary Operational Environmental Satellite) observations and described, for example, in Jackman et al. (2011). The energy range for protons is 1-300 MeV, and thus the direct atmospheric impact takes place at altitudes above ≈ 10 hPa (e,g, Turunen et al., 2009, Fig. 3).
Proton fluxes were only included for energies up to 300 MeV. However, Jackman et al. (2011) show that higher energies contribute very little to the total ozone or NO y impact. We only consider protons, and not the X-ray flares that are associated with some SPEs -in the polar regions we can assume that the proton effect is the dominant one, at least for large SPEs (Pettit et al., 2018).

Analysis methods
For our statistical analysis, SPEs were selected using proton flux data from satellite-based GOES observations (available from https://www.ngdc.noaa.gov/stp/ satellite/goes/index.html, last access: 9 July 2020). An event was selected if the peak proton flux exceeded 100 particle flux units (pfu), with pfu defined as the 5 min average flux in units of particles cm −2 s −1 sr −1 for protons with energy larger than 10 MeV. To avoid obscuring the signal in composite means by introducing duplicates of major events in days preceding the zero epoch date, any event following a larger event onset within 7 d was excluded from the analysis. In total, 11 such events were removed. Following these criteria, 66 events given in Table A1 were identified in the simulation period. The seasonal distribution of the selected SPEs is shown in Fig. 1. As SPEs are sporadic, the distribution is somewhat uneven, with the Northern Hemisphere winter months (December-February) underrepresented. Of the 8 largest SPEs (above 10 000 pfu) in the analysis, five occurred in October or November, compared to only one each in March, July and September.
For each of the 66 events, zero epoch day (D 0 ) was defined as the first day when the proton flux exceeded 10 pfu. A 90 d epoch period around the zero epoch ([D 0 −29, D 0 +60]) was selected from the daily mean time series for the following constituents: O 3 , HO x (OH + HO 2 ), Cl x (Cl + ClO), HNO 3 , NO x (NO + NO 2 ) and H 2 O. Definitions of abbreviations used in the following text for the statistical quantities from REF and WD simulations are shown in Table 1.
To assess the anomalies REF and WD associated with the SPEs, daily climatological mean mixing ratios REF and WD were deducted from respective epoch mixing ratios REF and WD. Figure 2 shows the composite means of WD for the northern polar cap (60-90 • N, area-weighted). Note that due to the seasonal distribution of the events some seasonal sig- nals remain visible in the composite means, e.g., a slight increase in mesospheric O 3 during the 90 d epoch period. It can also be seen that for some constituents, especially HNO 3 , mixing ratio enhancements following the events are large enough to visibly affect the daily climatology.
For each constituent, we calculate the composite mean of difference µ d between the epoch time series and the corresponding climatological time series  Fig. 1). Variance of the composite means was approximated by bootstrapping, i.e., by resampling the random selection of dates 1000 times with replacement, with N = 66 for each sample. Bootstrapping was chosen for variance estimation due to its independence from the shape of the distribution; i.e., it does not require a normal distribution. The standard deviation of the variance was then used to identify robust responses, i.e., SPE-driven changes that are clearly larger than the normal variability.
The SPE-driven anomalies are expected to occur following the SPE onset (zero epoch). However, in Sect. 3 we will show that the results sometimes have also anomalies that extend over the full epoch period. These anomalies are in most cases related to solar cycle variability, due to the SPE sampling over-representing the solar cycle maximum years when compared to climatology. Composite mean SPE ionization for the analyzed epochs is shown in Fig. 3. Daily composite mean peak ionization is 891 ion pairs cm −3 s −1 , which is roughly similar in magnitude to peak ionization associated with January 2005 SPE. The largest individual events in the time series feature peak ionizations about 10 times higher, for example peak ionization of 8534 ion pairs cm −3 s −1 for the "Halloween" SPE. In addition to the main ionization peak within ca. 5 d following the onset of the SPE, secondary ionization peaks are evident in the epoch mean, due to closely separated SPEs. While these secondary peaks are much weaker (by roughly an order of magnitude) than the main peak, they can still cause noticeable anomalies in some of the analyzed constituents.
In order to evaluate the contribution of the improved ion chemistry of WACCM-D, the same epoch analysis was made for both the REF and the WD simulations. The anomalies from the WACCM-D simulation ( WD) were then compared to anomalies from the standard WACCM ( REF) (Figs. 4 and  5). We considered the difference of composite means of differences in mixing ratio between the two simulations for the   Statistically robust anomalies, i.e., those larger than 2 times the standard deviation of the bootstrap variance, are shown in color and discussed below.
Several constituents show the most pronounced anomaly immediately following the event onset (zero epoch), which in the figures is marked with a black vertical line. For some constituents, most notably H 2 O and Cl x , persistent anomalies extending throughout the epoch period are seen as well. These anomalies can be considered to represent the solar cycle signal because the distribution of SPEs concentrates on the years around the solar maximum. Some anomalies are clearly affected by the response to closely separated SPEs, i.e., when the separation is smaller than 60 d. These are most notable in NO x and HNO 3 in NH around days −15, +22 and +45 with respect to zero epoch, corresponding to the secondary peaks of SPE ionization (see Fig. 3).
Ozone (Figs. 4 and 5, panel a) decreases between 1 and 0.01 hPa immediately following the SPE, lasting about 5-10 d. This is caused by catalytic destruction from enhanced HO x . Ozone loss also takes place around 1 hPa, driven by NO x and Cl x . The depletion around 1 hPa lasts much longer, driven by the NO x increase, with maximum ozone decrease seen about 30 d after the event onset. In contrast, an increase in ozone is seen throughout the period near the secondary ozone maximum above 0.01 hPa, more robustly in the Southern Hemisphere. This increase is linked to enhanced atomic oxygen production by O 2 photolysis in solar maximum conditions (see also . In the stratosphere below 5 hPa the ozone response is not consistently robust. However, an intermittent increase, likely caused by solar cycle effects, is seen throughout the epoch period in this region. Thus the decrease in the total ozone column, as reported by Denton et al. (2018), is not seen in our analysis. It should be noted that the highest energy protons (E > 300 MeV), which can directly impact the lower stratosphere, are not included in the particle forcing used in our simulations. However, it is not clear whether inclusion of these protons would change the response because their fluxes should be too low to cause a significant extra response .
A strong, short-lived increase in HO x is seen below 0.01 hPa (Figs. 4 and 5, panel b) immediately following the SPE, reaching down to the 10 hPa level, with the strongest and the most robust response in the mesosphere. Due to the short chemical lifetime of HO x , there is no longer-term anomaly beyond 10 d from the SPE onset. The later peaks just below 0.01 hPa are caused by SPEs occurring within 60 d of zero epoch of analyzed SPE. Above 0.01 hPa, a small negative anomaly is present throughout the period, related to solar maximum through less water vapor in the region and lowering of the HO x peak altitude. The persistent positive anomaly at just below 0.01 hPa, seen in the SH, is consistent with lowering of the HO x peak altitude. Factors affecting HO x peak altitude during solar maxima are the decrease in H 2 O and the increase in Lyman-α photolysis, balancing at the level of mesopause. Cl x mixing ratios (Figs. 4 and 5,panel c) are enhanced between 1 and 0.01 hPa, these SPE-driven increases last up to about a week. Below 1 hPa, at the altitude of Cl x mixing ratio maximum, Cl x is a reduced. This decrease is roughly consistent and continuous throughout the epoch period and so is likely, at least in part, connected to the solar cycle rather than individual events. However, the magnitude of the negative anomaly increases following the onset of the SPE, with maximum value in the SH reached after 30 d. A short-lived decrease in the NH polar upper stratosphere ClO during the January 2005 SPE has previously been reported both in satellite observations and models (Damiani et al., 2012;Andersson et al., 2016). In the short term decrease is likely caused by conversion of ClO to HOCl in the presence of enhanced HO 2 , while in the longer term, the decrease is probably connected to the descent of NO x and subsequent formation of ClONO 2 . Here we look at the polar average, but it should be noted that in the presence of SPE-generated HO x , the response of Cl x in the polar atmosphere depends strongly on the sunlight conditions and varies across latitudes and altitudes . For example, in the polar night stratosphere ClO is converted to HOCl in reaction with HO 2 , while in the daytime HOCl is converted to ClO by OH.
A large increase in HNO 3 (Figs. 4 and 5, panel d) is seen up to 10 d after the SPE onset, with a maximum between 1 and 0.01 hPa. The absence of a longer-term effect is due to the fast photodissociation of HNO 3 in sunlit conditions. Similar to HO x , the later peaks around days 20 and 45 in the NH are caused by nearly coincident SPEs. A longer-term enhancement persists for 20-30 d after SPEs below 1 hPa.  Fig. 4). The x axis shows the number of days before and after the event (day 0, solid black line) and the y axis pressure levels in the model (hPa, a, d) and approximate altitude (km, c, f). Note that the color scale for NO x panel is logarithmic, and all contours shown in that panel indicate positive difference.
The anomalies are weaker in the SH than in the NH, due to the strongest SPEs mostly occurring during NH winter (see Fig. 1) when there is less solar radiation and HNO 3 photodissociation.
SPE-driven enhancements of NO x (Figs. 4 and 5, panel e) reach down to 1 hPa level and below. These enhancements start to diminish after about 10 d, but robust longerterm enhancements can be seen for several weeks afterwards. This is especially clear at lower altitudes (around 1 hPa), where transport from above contributes to persistence of enhancements. The response is different in the Northern and Southern Hemispheres due to the seasonal distribution of the events. In the NH, several strong, closely separated wintertime events show up as recurring anomalies, which are then transported to lower altitudes. The SH response is less variable, with small but robust enhancements seen through the two months following the SPE onset. In addition to the SPEdriven anomaly, there is a persistent, underlying anomaly seen above about 10 hPa from enhanced NO x production in the lower thermosphere by solar EUV and EPP during solar maximum times (see also . During an SPE, H 2 O would be expected to decrease when it is converted to HO x in ionic reactions. In our analysis, H 2 O (Figs. 4 and 5, panel f) shows a consistent and statistically robust negative anomaly throughout the period ana-lyzed, except in the SH lower stratosphere. There is some indication of a stronger anomaly following the SPEs above 0.1 hPa, i.e., at altitudes where HO x increases. But in general the anomaly is mostly related to the solar cycle through the increased Lyman-α photodissociation during solar maximum decreasing H 2 O in the mesosphere-lower thermosphere (see also Schmidt et al., 2006;. Figures 6 and 7 show the difference between epoch responses in the NH and SH, respectively, from WACCM-D and the standard WACCM. As the specified dynamics of the two simulations are identical, differences seen in these figures arise from the addition of the D-region ion chemistry. Note that the robustness threshold for these figures is set to 1 standard deviation rather than the 2 standard deviations used in Figs Fig. 5). The x axis shows the number of days before and after the event (day 0, solid black line) and the y axis pressure levels in the model (hPa, a, d) and approximate altitude (km, c, f). Note that the color scale for NO x panel is logarithmic, and all contours shown in that panel indicate positive difference.

Effects of D-region ion chemistry
Ozone (Figs. 6 and 7, panel a) depletion is larger in WACCM-D for about 5 d after the onset at around 0.5-1 hPa. This effect is very similar in both hemispheres. Both the duration and altitude range of this extra ozone loss clearly correspond to enhanced Cl x , while there is no clear connection to NO x or HO x differences. Less ozone depletion takes place in WACCM-D around the 0.01 hPa level, connected to the clearly lower HOx enhancement in WACCM-D. An anomaly of similar magnitude is observed in both hemispheres, although it is only robust in SH. The reduced HO x response in WACCM-D is a result of the lower water vapor amount, and thus lower HO x production, than what is assumed when HO x production is parameterized . Cl x is enhanced in both hemispheres above 1 hPa in WACCM-D, while there is only a relatively weak increase in the standard WACCM. The reasons for the enhancement are the ion reactions included in WACCM-D that convert more HCl to active chlorine species Cl and ClO (Winkler et al., 2009) than the neutral gas-phase reactions which are included in both standard WACCM and WACCM-D. Enhancements relating to secondary ionization peaks can also be seen, most clearly around day 40. The negative anomaly of Cl x below 1 hPa seen in Figs. 4 and 5 is not present here, again implying that this anomaly is not predominantly due to ion chemistry.
The HNO 3 enhancement is clear and strong in WACCM-D, as ion chemistry produces it from other NO y species. Standard WACCM, like other models using an EPP lookup table parameterization, is known to underestimate the HNO 3 mixing ratios when compared to observations (Jackman et al., 2008;Funke et al., 2011;Andersson et al., 2016). A discrepancy between WACCM-D and standard WACCM HNO 3 enhancements is clearly seen in Figs. 8 and 9. Due to the dramatic increase in HNO 3 following SPEs, peaks from other, close-by events are also clearly seen in Fig. 8.
In the mesosphere, the NO x enhancement is much larger in WACCM-D than in standard WACCM. Below 0.1 hPa there are some signs of an increase in transported NO x , although this effect is not statistically robust. Increased NO x is consistent between the hemispheres. Andersson et al. (2016) attribute the increased NO x to enhanced production above 80 km (ca. 0.01 hPa), mostly due to the reaction O + 2 + N 2 → NO + + NO, and consequent transport to lower altitudes. Differences in H 2 O anomalies (Figs. 6-9, panel f) between WACCM-D and the standard WACCM are seen above 1.0 hPa but they are almost entirely overcome by the variance, which indicates that negative anomalies in Figs. 4 and 5 are not primarily connected to ion chemistry reactions but are a solar cycle signal. However, a short-lived negative difference, robust to a 1σ level, can be seen at 0.01 hPa immedi-8932 N. Kalakoski: SPE impact on atmospheric composition ately following the onset of SPE, even though the co-located HO x production in WACCM-D is smaller. This negative difference is more clearly seen in Figs. 8 and 9, where we see that mesospheric H 2 O anomalies are very consistent between the simulations, except within 10-15 d following the onset of SPE, where the WACCM-D negative anomalies are larger.

Effect from individual events
The SPE effects are dependent on the background atmosphere and season of the year. Here we discuss this by looking at the effects of individual events at selected altitudes. This is also an additional check on the robustness of our analysis. Figure 10 shows the individual time series of differences of epochs from the climatology around 0.05 hPa pressure level (calculated as the mean of three pressure levels centered at 0.04 hPa). SPEs used in this figure are screened for closely separated events, i.e., events which are followed by a larger event within 7 d are not shown. Individual events are sorted from top to bottom by the maximum observed proton flux, with strongest events on top of the figure. Dashed horizon-tal lines show the cut-offs for 1000 pfu (top line) and 100 pfu (bottom line) SPEs.
As expected, the largest SPEs cause a stronger response in O 3 , HO x and NO x . The response to events is more pronounced in the NH due to the seasonal distribution of SPEs; i.e., there are more events in NH winter. This is especially apparent for NO x , as its lifetime is heavily influenced by the amount of sunlight available. On visual inspection, the response is dominated by the largest SPEs, with anomalies following the onset of SPEs being visually indistinguishable for the weakest events. Below the 100 pfu threshold, any anomalies following the event onset are dominated by larger events occurring within the analysis period, rather than by the analyzed event itself. Thus, the decision to only analyze larger events seems reasonable because the events smaller than 100 pfu would only add a substantial amount of noise to the analysis.

Conclusions
We present an analysis of the chemical impacts of SPE on the middle atmosphere using simulations from WACCM- The aim of our analysis is to study the impact and improvement from a detailed ion chemistry scheme, which is done by comparing WACCM-D results to those from the standard WACCM scheme which used a simple parameterization of HO x and NO x production. Instead of analyzing individual SPE events, we present a statistical analysis of the 66 largest SPEs over the years 1989-2012 to allow for more general conclusions in a range of different background atmospheric conditions and with the observed seasonal distribution of SPEs. This statistical approach allowed for the incorporation of smaller events which would have been difficult to analyze separately. Analysis was performed by means of a superposed epoch method with bootstrapping to identify statistically robust responses.
Statistically robust SPE signals were present in WACCM-D for all analyzed species other than water vapor (O 3 , HO x , Cl x , HNO 3 , NO x ). These signals are qualitatively in line with previously published modeling results and satellite observations considering individual events. In addition, our analysis shows longer-term, solar-cycle-type signals for several species due to SPEs predominantly occurring during solar maximum years. In general, the responses are consistent between the two hemispheres.
Compared to the standard WACCM, WACCM-D provides a solution to some known shortcomings in chemistry-climate models, particularly in the case of HNO 3 and Cl x response to SPEs. Ozone loss following the SPEs was enhanced at 0.5 hPa in WACCM-D. As there is no corresponding increase in HO x or NO x at this altitude, this loss is connected to increased conversion of HCl to reactive Cl x species by the ion reactions (Winkler et al., 2009), now accounted for in WACCM-D. Conversely, less ozone loss took place around 0.01 hPa in the SH, due to less production of HO x , in agreement with Andersson et al. (2016). NO x production in the mesosphere is enhanced by the inclusion of ion chemistry and subsequent downward transport, extending the effect to lower altitudes for up to 20 d following the SPE onset. Ion chemistry was found to be crucial for characterization of HNO 3 production following the SPEs, which confirms the previous results from 1-D chemistry models (Verronen et al., 2008;Verronen et al., 2011). Some effect was also found in H 2 O in an altitude region near the mesopause (around 0.01 hPa or 85 km). This, however, has a large 11-year solar cycle signal, partially masking the SPE impact. In summary, incomplete representation of EPP-driven chemistry in simulations, for example regarding the HNO 3 , has been recognized as one of the outstanding questions in understanding the solar influence in the middle atmosphere and below (Jackman et al., 2008;Funke et al., 2011). Our results clearly show the importance of D-region ion chemistry in capturing the effects of energetic particle precipitation in the mesosphere-lower thermosphere simulations. Improved global modeling with ion chemistry, such as that in WACCM-D, provides an important tool for interpretation of a wider range of satellite-based observations of neutral species and allows for global studies of ionospheric D region.
Author contributions. NK, PTV and MES planned the research and simulation. MES carried out the simulations. NK carried out the analysis and led the writing of the paper. All authors contributed to the discussion and the writing of the final paper.
Competing interests. The authors declare that they have no conflict of interest.
Acknowledgements. This material is based upon work supported in part by the National Center for Atmospheric Research, which is a major facility sponsored by the National Science Foundation under cooperative agreement no. 1852977. Work was carried out as a part of International Space Science Institute (ISSI) project "Space Weather Induced Direct Ionisation Effects On The Ozone Layer".
Financial support. DRM has been supported by NSF (Award no. 1650918, "Collaborative Research: CEDAR -Quantifying the Impact of Radiation Belt Electron Precipitation on Atmospheric Reactive Nitrogen Oxides (NO x ) and Ozone (O 3 )"). Antti Kero is funded by the Tenure Track Project in Radio Science at Sodankylä Geophysical Observatory.
Review statement. This paper was edited by Gabriele Stiller and reviewed by three anonymous referees.