Will Climate Change Impact Polar NOx Produced by Energetic Particle Precipitation?

Energetic electron precipitation (EEP) is an important source of polar nitrogen oxides (NOx) in the upper atmosphere. During winter, mesospheric NOx has a long chemical lifetime and is transported to the stratosphere by the mean meridional circulation. Climate change is expected to accelerate this circulation and therefore increase polar mesospheric descent rates. We investigate the Southern Hemispheric polar NOx distribution during the 21st century under a variety of future scenarios using simulations of the Whole Atmosphere Community Climate Model (WACCM). We simulate stronger polar mesospheric descent in all future scenarios that increase the atmospheric radiative forcing. Polar NOx in the upper stratosphere is significantly enhanced in two future scenarios with the largest increase in radiative forcing. This indicates that the ozone depleting NOx cycle will become more important in the future, especially if stratospheric chlorine species decline. Thus, EEP‐related atmospheric effects may become more prominent in the future.


Introduction
Energetic particle precipitation (EPP) is an important source of ionization in the upper atmosphere (Brown, 1966;Sinnhuber et al., 2012). These particles consist of electrons and protons precipitating from the solar wind and the magnetosphere. The peak ionization altitude depends on the type of particle and its initial energy (Vampola & Gorney, 1983;Mironova et al., 2015). Energetic electron precipitation (EEP) from the magnetosphere consists of low-energy (auroral) electrons with ionization peaking in the lower thermosphere (Turunen et al., 2009), and medium energy electrons with peak ionization in the upper mesosphere (Nesse Tyssøy et al., 2016). Auroral electrons (<40 keV) come from the magnetospheric plasmasheet and are accelerated by substorm activity, whereas medium-energy electrons originate from the radiation belts and are accelerated by wave-particle interactions during geomagnetic storms (Horne et al., 2009). EEP is confined to the high geomagnetic latitudes, where it forms an oval shape pattern circling the magnetic pole.
Other important energetic particle sources are solar proton events (SPE) and galactic cosmic rays (GCR). SPE are sporadic events of high-energy proton injection directly from the solar wind. They can penetrate over the entire polar cap, down to the upper stratosphere (Jackman et al., 2008). GCR come from outside of the solar system and have enough energy to precipitate down to the troposphere (Usoskin & Kovaltsov, 2006). However, their total number and average ionization rate are several magnitudes smaller than that of EEP or SPE in the mesosphere and upper stratosphere (Mironova et al., 2015).

10.1029/2020GL087041
Geomagnetic disturbances have been directly observed on the surface since the nineteenth century (Svalgaard & Cliver, 2010) and can be used as a proxy for EEP. For example, globally averaged geomagnetic activity index Ap roughly corresponds to the ionospheric electric currents produced by auroral electrons (Østgaard et al., 2002). Recently, this has been used to reconstruct EEP activity since the midnineteenth century (Matthes et al., 2017).
EEP in the thermosphere, sometimes called auroral EEP since it can produce strong airglow, is an important source of nitrogen oxides (NOx=N+NO+NO 2 ) (Solomon et al., 1982;Smith-Johnsen et al., 2018). SPEs and medium energy (>40 keV) EEP produce NOx and hydrogen oxides (HOx=H+OH+HO 2 ) in the mesosphere (Solomon et al., 1981;Sinnhuber et al., 2012) (in case of SPEs also in the upper stratosphere (Jackman et al., 2008)). HOx is short-lived, and dissipates fairly quickly (Verronen et al., 2011). NOx, on the other hand, has a prolonged lifetime in the polar night mesosphere where it is not photolyzed. Consequently, substantial increases in polar NOx are observed in the winter time (Funke et al., 2014).
EEP-related polar NOx enhancements are also observed in the stratosphere, but with a time lag relative to the mesospheric enhancements (Randall et al., 2006). This is because the dominating mean atmospheric flow is a poleward and downward circulation in the high latitudinal mesosphere during winter (Andrews et al., 1987), and this circulation slowly transports mesospheric NOx into the stratosphere. This is often referred to as the indirect effect of EEP-NOx (Randall et al., 2006).
Both HOx and NOx are effective ozone destroyers, which leads to observable ozone depletion in the mesosphere following elevated geomagnetic activity, and with some lag in the stratosphere (Andersson et al., 2018). Ozone is a radiatively active gas, and thus EEP-related chemical changes can impact the dynamics of the middle atmosphere (Baumgaertner et al., 2011).
Both observations (Fu et al., 2015) and climate model predictions (Butchart et al., 2006) indicate that the mean meridional circulation is accelerating due to the climate change. Since NOx varies with the mean meridional circulation, it is reasonable to ask how EEP-related NOx will vary in the future. This is the topic of this paper. Baumgaertner et al. (2010) have shown that during high levels of geomagnetic activity and with the most drastic future emission scenario, average winter polar NOx levels will be higher. We will study future variability of the EEP-NOx over the whole 21st century under four different future scenarios, each having the same moderate variable solar activity scenario. In section 2 we describe the climate model data and the statistical methods utilized. In section 3 we present how polar NOx depends on different drivers in the historical model run. In section 4 we present changes in polar NOx in model runs of the different Shared Socioeconomic Pathways (SSPs). Finally, our conclusions are presented in section 5.

Data and Methods
We use output from different simulations of a free-running version of the chemistry-climate model WACCM6 (Whole Atmosphere Community Climate Model) within CESM2. The model structure is explained in detail by Marsh et al. (2013) and updated in Gettelman et al. (2019). Five different simulations are analyzed: Historical simulation (three ensemble members) covering the period 1850-2014 (CMIP6 DECK simulations), and four different future scenario (CMIP6 ScenarioMIP: SSP1, SSP2, SSP3, and SSP5) simulations with one ensemble member each covering the period 2015-2100 (O'Neill et al., 2016). Different SSPs cover a wide range of actions by society, and not just varying carbon dioxide concentrations. The resulting radiative forcing increase of the climate system by 2100 relative to preindustrial era is projected to be 5.0 W/m 2 in SSP1, 6.5 W/m 2 in SSP2, 7.2 W/m 2 in SSP3, and 8.7 W/m 2 in SSP5. One can obtain details of the different SSPs in Riahi et al. (2017).
We focus on monthly mean zonal mean data, mainly in the Southern Hemisphere winter (May to September). The model has latitudinal resolution of roughly 0.94 • (192 bins) and altitude range from the surface up to 140 km (in 70 levels). In this study we focus on modeled volume mixing ratios of NOx from around the mesopause to the surface (0.01 to 1,000 hPa), as well as the vertical component of the residual circulation, w * (Andrews et al., 1987). In addition, we also obtain directly from the model output indices for ENSO (El Niño-Southern Oscillation; Nino3.4: detrended sea surface temperature averaged between 5 • S-5 • N and 170 • -120 • W), QBO (Quasi-Biennial Oscillation; zonal mean zonal wind in 5 • S-5 • N at 15 and 32 hPa) and volcanic activity (zonal mean aerosol optical depth 550 nm in the 90 • S-90 • N stratosphere).

10.1029/2020GL087041
The Ap geomagnetic index and SPE (ion-pair production rate at 1 hPa between 66 and 82 • S) are prescribed following the recommendations of Coupled Model Intercomparison Project 6 (CMIP6). This provides estimates of the solar activity before the space era, and a future solar forcing scenario (1850-2100 in this case) Matthes et al. (2017). Model runs also include variable solar spectral irradiance and GCR, but they are not considered in this study. As stated above, GCR produce some NOx in the troposphere and lower stratosphere, but the effect is small (less than 1% increase from solar maximum to minimum) above 10 hPa Jackman et al. (2016).
Multiple linear regression (MLR) is used to fit time series of NOx at each altitude and latitude in the historical run using six regressors: Ap index, SPE, ENSO, QBO (15and 32 hPa), and volcanic activity. Two QBO terms are used to account for the fact that the QBO is not constant with height (There is roughly a quarter cycle shift between the QBO at 15 and at 32 hPa). The cross correlations between each of these regressors for each month is shown in Tables S1 to S5 in the supporting information. All of the cross correlations are fairly small (between−0.20 and 0.22) and thus we can assume that the multicollinearity is not a big issue. MLR is calculated using the Cochrane-Orcutt method (Cochrane & Orcutt, 1949). In this method, the residual term is obtained as an autoregressive AR(1) process instead of normally distributed white noise as in normal linear regression. Accordingly, where i represents each explaining variable and the residual term is t = t−1 + e t with being the lag-1 autocorrelation of the residual ande t is normally distributed white noise. This can be written as i , and can be solved iteratively by first calculating the normal regression for (1) and estimating from the residual time series t . A regression is then performed using (2) with the value of obtained in the first step. After that a new estimate for is calculated from the new residuals. This can be iterated until converges. The significance of regression parameters can be calculated using the Student's t test, which would not be correct with a residual term that was autocorrelated (note that this is an idealization since AR(n > 1) processes might also occur). More details of the procedure can be obtained in Maliniemi et al. (2018). The Ap index is averaged over two months (Apr/May, May/Jun, Jun/Jul, Jun/Jul, Jun/Jul), which occur before the NOx response (May, Jun, Jul, Aug, Sep) to obtain the indirect effect due to NOx descending from the thermosphere and upper mesosphere (Funke et al., 2014). Other regression terms have the same 1 month average as the NOx response. Regression results for the historical run are averages over the three ensemble members.
Statistical significance for the differences in NOx between 2090-2100 and 2015-2025 for the different SSPs are calculated applying a Monte Carlo method: we take two random 11-year time periods from 2015-2100 and calculate the difference in each latitude/height bin. This is performed 10,000 times and original value (difference of 2090-2100 and 2015-2025) is compared to the distribution of these 10,000 repetitions to obtain fraction of more extreme differences (both tails). This fraction then represents the p value in each bin.
Because we present results over several latitudes and pressure levels, we have a multiple hypothesis testing situation. Thus, simply presenting significance in each bin based on individual hypothesis testing would lead to overestimation of the true number of rejected null hypotheses. This is because of the interpretation of the p value, and the dependency of the neighboring grid points, that is, the spatial autocorrelation (Wilks, 2016). To overcome this issue, we have obtained the false discovery rate. It adjusts the used p value limit of individual hypothesis tests to take into account the spatial autocorrelation. After the procedure, the p value can be interpreted as being the probability to obtain a false rejection of null hypothesis. Details of the method can be found in Wilks (2016). P value limit in our case is 0.05, which then represents the rejection of global null hypothesis with 95% probability, when in fact local p value limits (in individual bins) are notably smaller.
The smoothly varying long-term variations shown in section 4 are calculated for upper stratospheric NOx, Ap index, and mesospheric vertical residual circulation. This is done using the LOWESS-method (LOcally WEighted Scatterplot Smoothing) applied with a 31-year window (Cleveland & Devlin, 1988).  Figure 1 presents the May to September Southern Hemisphere polar NOx response to 1 standard deviation increase in the Ap and SPE terms used in the MLR (Ap index time series normalized as[x − ]∕ and SPE asx∕ ). One can see that the SPE-related NOx peaks in the lower mesosphere but the distribution extends also to the upper mesosphere and upper stratosphere. The altitudinal distribution of changes in volume mixing ratio are roughly constant across the winter months. It is slightly larger in May (approximately 50%) and smaller in June, but this is probably due to the reasonably short time period, 1850-2014, which is not long enough to distribute large SPEs evenly over the different months.

Historical Polar NOx Response to Different External Forcings
The relative NOx response to energetic electron precipitation, represented by the Ap index, peaks in the upper mesosphere in the early winter but descends to the upper stratosphere later in winter. Again, this is consistent with the fact that NOx has a long chemical lifetime and is transported from the thermosphere/upper mesosphere during winter by the residual circulation. The NOx enhancement to 1 standard deviation increase of Ap index is roughly 20%. Both the Ap index and SPE responses are significant over large-altitude range in all months.
SPEs and periods of elevated geomagnetic activity (high Ap) produce significant NOx increase in the high latitudes and in the mesosphere/upper stratosphere. One can see the absolute NOx response to each regressor over the whole Southern Hemisphere between May and September in Figure S1 in the supporting information. Other explaining variables produce only small and insignificant response to NOx in the high latitudinal mesosphere and upper stratosphere. Volcanic activity produces a significant negative response mainly in the lower stratosphere, which partly extends to high latitudes. The response to ENSO is mainly negative in the tropical troposphere. The two QBO terms produce responses in the tropical and subtropical  Figure 2. In the historical period large spikes in the NOx amount are directly related to SPE activity, where as the overall (decadal) level of NOx follows the Ap index quite well. Polar mesospheric w * is relatively constant over the historical run. This indicates that the majority of the decadal polar NOx variability in the historical period is related to EEP, rather than atmospheric circulation.

Future Polar NOx Ditributions in Different SSPs
Figure 2 also presents stratospheric polar NOx and w * in simulations under different SSPs during 2015-+2100. There is a similar and slightly positive trend in NOx in all SSPs up to the 2040s, at which point all SSPs show a substantial increase due to the several SPEs occuring in the 2050s. Note that SPE activity after 2014 (and also prior 1963) is based on satellite measurements of proton fluxes during 1963-2014, which is then (quasi randomly) projected according to the overall solar activity level (Matthes et al., 2017). Thus, the peak in the SPE frequency around 2050s is really a feature of the forcing scenario. The different SSPs differ more toward the latter half of the 21st century. NOx in SSP1 returns roughly to the pre 2050s level, whereas in the other SSPs there is more upper stratospheric NOx than in historical run or before the 2050s. This increase in SSP2, SSP3, and SSP5 cannot be explained by the occurence of SPEs (absence of major events during 2060-2100) or by the Ap index (which is at the same level or lower than during the twentieth century).
However, the vertical descent in the polar mesosphere shows notable differences across the varying SSPs. All SSPs produce an increasing trend in the descent rate between 2015 and 2050, after which SSP1 and SSP2 starts to level out. In SSP3 and SSP5 the trend continues to the end of the century, especially in SSP5 where Geophysical Research Letters 10.1029/2020GL087041 the descent rate does not appear to level off. The trends in NOx are consistent with the w * trends, if one excludes the period of high SPE activity around 2050. Figure 3 presents the monthly absolute differences in NOx between averages over the time periods 2090 to 2100 and 2015 to 2025 for SSP5. One can see that the largest differences are obtained in the polar mesosphere and upper stratosphere during winter, especially in the Southern Hemisphere. NOx differences are greatest over the southern pole in all winter months, which indicates that it is related to the increased descent in the polar mesosphere.
In the north there is also an increase in the polar NOx, but it is restricted mostly to the upper mesosphere. The polar atmosphere in the north experiences much more variability from one winter to the next due to the more planetary wave activity (Salby & Callaghan, 2003), and the occurrence of sudden stratospheric warmings (SSWs) (Manney et al., 2005). Related to this, we also checked if SSWs were obtained in any of the model runs. None were simulated in the Southern Hemisphere in the historical or any of the SSP runs, whereas in the north an average of 3-6 events per decade were obtained in all runs. Here we used the SSW definition of a reversal of daily zonal mean zonal wind at 60 • latitude and 10 hPa during three central winter months. Thus, Northern Hemisphere NOx responses below the upper mesosphere is likely affected by SSWs and accompanying elevated stratopause events, as previously shown by Holt et al. (2013).
One can see the global distribution of NOx differences between 2090-2100 and 2015-2025 in all SSPs during southern winter in Figure S2 in the supporting information. Increased polar NOx in SSP5 is also seen SSP3, whereas in SSP1 and SSP2 there are hardly any significant differences in the polar region. There are also other trends present in the atmosphere in different SSPs. Tropospheric NOx has negative trends in all SSPs, and equatorial stratospheric NOx has positive trends in all SSPs (it is strongest in SSP3). However, these are substantially smaller, and presumably not the cause of polar NOx trends because the dominating flow is mostly downward in the winter polar mesosphere and upper stratosphere.
Finally, the relative differences in polar cap NOx between 2090-2100 and 2015-2025 for all SSPs are shown in Figure 4. One can see that in all months the significant NOx increase is largest in SSP5, followed by SSP3, SSP2, and SSP1, though mostly insignificant in SSP2 and SSP1. This supports the NOx trends in Figure 2. Looking at the distributions in Figure 4, SSP3 and SSP5 match quite well the Ap index induced distributions in Figure 1. This supports the interpretation that the increase of NOx originates from geomagnetic activity/energetic electron precipitation above the stratopause that is brought down by enhanced vertical descent.

Conclusions
In this letter we have shown that winter polar NOx is strongly related to the geomagnetic activity and solar proton events. Polar NOx response to SPEs is relatively constant over different winter months peaking in the lower mesosphere. Polar NOx responses to Ap index variability changes over the winter and spring, peaking in the upper mesosphere in the early winter and descending to the upper stratosphere in spring.
Polar NOx is significantly increased in SSP3 and SSP5 toward the end of the 21st century. This increase cannot be explained by the occurrence of strong SPEs (which are mostly absent between 2060 and 2100) or the level of geomagnetic activity (which is not higher than during the twentieth century) in these scenarios. However, the polar mesospheric descent rate is greatly increased in SSP3 and SSP5, whereas in SSP1 and SSP2 w * trends level off toward the end of the 21st century.
Distribution of increases in polar NOx in SSP3 and SSP5 (Figure 4) is very similar to the distribution of polar NOx related to geomagnetic activity determined by MLR in the historical run (Figure 1). This strongly indicates that the future polar NOx increase in the stratosphere are a consequence of a stronger EEP indirect effect that originates in the thermosphere and upper mesosphere.
Our results for SSP5 are similar to those obtained by Baumgaertner et al. (2010) who looked at the response of NOx under a single future scenario (SRES A2, which has a near doubling of CO 2 in 2100) and specification of the elevated geomagnetic forcing taken from 2003. Our results suggest that the future increase of polar NOx do not necessarily need high level of geomagnetic activity, but is a direct consequence of circulation changes. It also shows that increases in polar NOx can also be obtained in the more moderate future scenarios (e.g., SSP3).
These results give interesting implications to the projections of the future atmospheric state. Polar NOx enhancements by EEP can also impact atmospheric dynamics (Salminen et al., 2019) due to the ozone depletion by catalytic NOx cycle (Andersson et al., 2018). CFC-emissions have been successfully reduced following adoption of the Montreal Protocol (Velders et al., 2007) and the Antarctic ozone hole is already showing signs of recovery (Solomon et al., 2016). Our results imply that the catalytic NOx cycle will become relatively more important in the future, while the ClOx cycle is diminishing. Thus, the EEP-related atmospheric effects will become stronger when the climate change proceeds, regardless of the level of EEP activity.