Influence of Enhanced Planetary Wave Activity on the Polar Vortex Enhancement Related to Energetic Electron Precipitation

The Northern polar vortex experiences considerable interannual variability, which is also reflected to tropospheric weather. Recent research has established a link between polar vortex variations and energetic electron precipitation (EEP) from the near‐Earth space into the polar atmosphere, which is mediated by EEP‐induced chemical changes causing ozone loss in the mesosphere and stratosphere. The most dramatic changes in the polar vortex are due to sudden stratospheric warmings (SSWs). Enhanced planetary wave convergence and meridional circulation may cause an SSW, a temporary breakdown of the polar vortex. Here we study the relation of SSWs to the atmospheric response to EEP in 1957–2017 using combined ERA‐40 and ERA‐Interim reanalysis data and geomagnetic activity as a proxy of EEP. We find that the EEP‐related enhancement of the polar vortex and other associated dynamical responses are seen only during winters when an SSW occurs and that the EEP‐related changes are observed systematically slightly before the SSW onset.


Introduction
During Arctic winter, the polar stratosphere remains mostly in darkness and cools radiatively. The resulting difference in air temperature between the pole and lower latitudes leads to the formation of a strong westerly zonal wind, the polar vortex (e.g., Andrews, 2010). A polar vortex is also formed in the Southern Hemisphere during Antarctic winter. The strength of the polar vortex especially in the Northern Hemisphere varies considerably during the winter and from one year to the next. Variations in the polar stratospheric temperature and, thus, the strength of the polar vortex are correlated well with the variations of planetary waves propagating from the troposphere to the stratosphere (e.g., Salby & Callaghan, 2002). Converging wave activity imparts enhanced drag on the westerly flow, enhances mixing of air between polar and lower latitudes, and Polar stratospheric temperatures and winds also depend on the amount of ozone. Ozone is a radiatively active trace gas, which can absorb solar UV radiation and thermal long-wave radiation from below and also emit long-wave thermal radiation (Brasseur & Solomon, 2005). Observations and modeling studies (e.g., Baumgaertner et al., 2011;Langematz et al., 2003;Sinnhuber et al., 2018) have shown that ozone loss during winter leads to net radiative warming (less long-wave radiation emitted to space) in the polar upper stratosphere/mesosphere and to a net cooling in the lower stratosphere (less long-wave radiation absorbed from below). These ozone loss-associated temperature changes in the stratosphere are also expected to enhance the polar vortex by thermal wind shear balance.
Energetic particle precipitation from the Earth's magnetosphere contributes to ozone loss (Callis et al., 1991;Randall et al., 2005;Rozanov et al., 2005;Turunen et al., 2009) by creating NO x /NO y and HO x family of chemical species in the mesosphere, which then catalyze ozone loss (Grenfell et al., 2006;Solomon, 1999), either directly in the mesosphere or indirectly somewhat later, after they are transported down to the stratosphere during winter (Funke et al., 2014a). The descent of NOx/NOy occurs every winter, and significant amounts are seen in the stratosphere at 1-to 10-hPa altitudes typically already from December onward . It has been shown that energetic electron precipitation (EEP) causes significant ozone loss by 5-20% in the mesosphere and stratosphere during wintertime (Andersson et al., 2014(Andersson et al., , 2018. Ozone loss influences the polar stratospheric temperature and the polar vortex, and the changes in the polar vortex contribute to the wintertime weather variations on ground Douville, 2009;Kidston et al., 2015;Scaife et al., 2005).
Long-term observations indicate that interannual variations in wintertime weather related to EEP are modulated by the QBO. For example, Maliniemi et al. (2013) and Maliniemi et al. (2016) showed that the EEP-related wintertime climate variations on ground are clearly stronger and more significant in the easterly phase of QBO. Similar results for the stratospheric polar vortex were recently reported by Salminen et al. (2019) who showed that EEP-related enhancement of polar vortex is much stronger in the easterly QBO phase. They suggested that the observed QBO modulation might be related to the fact that easterly QBO is associated with stronger planetary wave activity and stronger meridional circulation/downwelling in the polar stratosphere, which might make the polar vortex more susceptible to wave-mean-flow interactions.
In this paper we study the effect of EEP on the stratosphere and consider in more detail how the changes associated with SSWs influence the way EEP affects the polar vortex dynamics during the winter. Since SSWs are associated with enhanced planetary wave activity and changes of meridional circulation, we will first briefly review in section 2 their overall connection to zonal wind and temperature. In this study we will utilize a composite climate data set constructed from ERA-40 and ERA-Interim reanalysis data, which together cover the time from 1957 to present. As a proxy of EEP for this time period, we shall use geomagnetic activity characterized by the Ap index. These data sets will be discussed in section 3, and their analysis methods will be discussed in section 4. We then discuss the Ap-related responses in the atmosphere and highlight the importance of SSWs for these responses in sections 5 and 6, respectively. The last section summarizes our results with concluding discussion and remarks.

Theoretical Considerations
In this work we consider the effects of planetary wave activity and meridional residual circulation on atmospheric dynamics. Planetary waves are characterized by the Eliassen-Palm (EP) flux vector. In a coordinate system, where the vertical coordinate is the air pressure p, the meridional and vertical components of the EP-flux vector in quasi-geostropic approximation are defined as where is latitude, r 0 = 6, 371 km is the Earth's radius, = 2Ω sin is the Coriolis parameter, = 2 ∕86, 400 s is Earth's angular velocity, v is meridional wind, u is the zonal wind, and is the potential temperature (Andrews et al., 1987;Edmon et al., 1980). Potential temperature is defined as = T(p 0 ∕p) R∕c p ≈ T(1, 000 hPa∕p) 0.286 . The primed quantities v ′ , u ′ , and ′ refer to differences of the quantities from their zonal averages. The overbar indicates zonal average. The divergence of EP-flux in a spherical coordinate system with p as the vertical component can be written as (3) Andrews et al. (1987) and Edmon et al. (1980) showed that one can write the meridional and vertical components of residual circulation asv * =v − p where w is the vertical velocity in units of Pa/s. Note that in these units, upwelling corresponds to negative, and downwelling to positive vertical velocity. Residual circulation is connected to the evolution of zonal-mean zonal windū by the zonal momentum equation. In the absence of small-scale friction, the quasi-geostropic zonal momentum equation can be expressed in transformed Eulerian mean formulation as (Edmon et al., 1980) This equation indicates that the acceleration/deceleration of zonal mean zonal wind is due to the zonal component of the Coriolis acceleration (first term on the right-hand side, which is related to meridional component of residual circulation) and convergence/divergence of planetary waves (second term on the right-hand side). One can also formulate in transformed Eulerian mean framework the thermodynamic equation describing the evolution of potential temperature as where Q represents diabatic heating rate, for example, due to absorbtion of solar radiation. Note that when pressure is used as the vertical axis, the potential temperature gradient ∕ p is typically negative. Thus, equation (7) indicates that in addition to Q, the potential temperature may increase (decrease) if downwelling increases (decreases). Equations (6) and (7) imply that the momentary values of zonal wind and potential temperature are obtained as time integrals of the various terms on the right-hand sides of the equations. For example, when considering the effect of planetary waves on the zonal wind at a particular moment of time, it is appropriate to compute the time integral of the EP-flux divergence. Past studies (e.g., Newman et al., 2001;Polvani & Waugh, 2004) have indeed shown that stratospheric zonal wind depends not on the momentary EP fluxes but rather on values integrated over a period of time corresponding roughly to the thermal relaxation timescale in the stratosphere, which is roughly 10-60 days depending on the altitude (e.g., Mlynczak et al., 1999).

Atmospheric Reanalyses Data
In this study we use the ERA-40 (Uppala et al., 2005) and ERA-Interim (Dee et al., 2011) atmospheric reanalysis data obtained from https://apps.ecmwf.int/datasets/. ERA-40 provides data from September 1957 until August 2002, while the ERA-Interim data used here span from January 1979 to March 2017. The data are offered in latitude-longitude grid of 2.5 • × 2.5 • resolution with 23 pressure levels in ERA-40 and 37 in ERA-Interim. The pressure levels are roughly logarithmically spaced between 1,000 hPa (surface) and 1 hPa (roughly the stratopause level). In this work we used 6-hourly data from both reanalyses to first calculate daily means of zonally averaged zonal wind and temperature. We also computed in 6-hr time resolution the components of the EP flux and its divergence (equations (1-3)), which are used for diagnosing planetary wave activity as well as the two components of the meridional residual circulation (equations 4 and 5). We also calculated the time-integrated daily values of EP-flux divergence term in equation (6) and the adiabatic heating rate (second term on the right-hand side of equation 7) as a sliding sum using a [−30, … , 0 day] window, which roughly corresponds to the thermal relaxation time in the middle stratosphere (Mlynczak et al., 1999). In other words, the integrated EP-flux divergence (and similarly the integrated adiabatic heating rate) at time t was calculated as a sum where t = 86, 400 s. Note, however, that the results of this paper considering these integrated quantities are not significantly dependent on the particular choice the integration window length, as long as it is on the order of magnitude of thermal relaxation time.
In this work we combine the ERA-40 and ERA-Interim into a composite data series. It should be noted that, although the data of the two reanalyses are fairly close to each other, there are some significant differences, for example, in temperature and the EP fluxes in various regions during the overlapping time period from January 1979 to August 2002. These differences arise from the differences in the assimilated data and the construction of the underlying model and data assimilation methodology (Dee et al., 2011;Uppala et al., 2005). To reduce the differences between the two data sets, we scaled the ERA-40 data to ERA-Interim level with the general methodology developed by Asikainen (2019). The scaling was done separately for each parameter and for daily resolution data. The monthly averages were then calculated from the composite daily data set.
We will briefly describe the main steps of the method here. First, the data sets are standardized after which the dimensions of both data sets are reduced by computing the principal components (PCs) and corresponding empirical orthogonal functions (EOFs) using the data from the common time period 1 January 1979 to 31 August 2002. Instead of all PCs, only a subset of the most important PCs is retained for further analysis, and the rest of the variability in the data is considered statistical noise, which is uncorrelated between the two data sets. The optimal number of PCs for both sets is determined by the iterative DINEOF (Data Interpolating Empirical Orthogonal Function) method (Beckers & Rixen, 2003). DINEOF computes PCs by setting aside a randomly selected portion of the data (in our case 10%) to be used for optimization. Optimal number of PCs is determined by estimating the values of the optimization data with successively increasing number of PCs and stopping at the point, where the mean-squared error does not significantly decrease anymore (relative change of the error is less than 1%). After finding the optimal number of PCs, the ERA-40 set of PCs is mapped to ERA-Interim level using canonical correlation analysis (CCA). The CCA finds sets of linear combinations of the PCs (canonical variables), which are maximally correlated between the data sets so that each canonical variable is uncorrelated with all other canonical variables. In other words, CCA finds matrices A and B so that where P ERA-Interim and P ERA-40 are matrices that contain the time series of the optimal number of PCs from ERA-Interim and ERA-40, respectively, and is a matrix of residual noise. An estimate of the ERA-Interim PCs can then be obtained as With equation (10), we can estimate what the ERA-Interim PCs would be in the entire ERA-40 time interval by using only ERA-40 data. This is done by projecting the ERA-40 data onto the EOFs based on the data only from the common time period. The estimated ERA-Interim-level PCs (equation 10) are then multiplied by the corresponding ERA-Interim EOFs, and the result is scaled and centered to match the ERA-Interim standard deviation and mean of the common time period (from which the EOFs and the mapping in equation 10 were determined). Finally, the remaining residuals in the ERA-40 data set are rescaled and centered to have the same standard deviation and mean as the residuals in ERA-Interim in the common time period. The composite data set is then formed by ERA-40 data from September 1957 to 31 December 1978 and ERA-Interim data from 1 January 1979 to 31 March 2017. Figure 1 shows an example of how the ERA-40 scaling performs in the upper polar stratosphere at 70 • N, 5-hPa height. One can see that in this region, the two reanalyses differ from each other significantly and even the annual variation in temperature has a significantly larger amplitude in ERA-40. The average difference between daily unscaled ERA-40 and ERA-Interim temperature in this region is −2.5 K, while the scaling removes this average difference entirely. The standard deviations of the ERA-40/ERA-Interim differences in the considered region are 3.2 K and 1.4 K, respectively. Thus, the scaling reduces the residual differences by more than a factor of 2. Overall, the scaling brings ERA-40 data very close to the ERA-Interim data and significantly improves the homogeneity of the composite. The polar upper stratosphere was selected as an example in Figure 1 because the differences between the two reanalyses appear to be largest there. In the lower polar stratosphere (e.g., at 50-hPa height), the reanalyses agree much better, and the scaling does not change the data as much (not shown).

Geomagnetic Activity
As a proxy of EEP into the atmosphere, we use geomagnetic Ap index obtained from the GFZ Helmholz-Zentrum Potsdam, Germany (https://www.gfz-potsdam.de/en/kp-index/). Ap index measures the range of geomagnetic variability in 3-hr intervals and is an average over 13 midlatitude geomagnetic observatories. The Ap index responds to many different electric current systems in near-Earth space, but the largest contributions to it come from the field-aligned currents and their closure currents in the ionosphere. These currents have largest disturbances during auroral substorms and geomagnetic storms, when also energetic particle acceleration and precipitation are maximized. In longer timescales geomagnetic activity and the associated energetic particle precipitation are dominantly driven by high-speed solar wind streams emanating from solar coronal holes (Asikainen & Ruopsa, 2016;Richardson & Cane, 2012). While their evolution is related to the solar cycle, they do not follow the sunspot number. Geomagnetic activity has local minima in sunspot minima and but maxima typically a few years after the sunspot maximum (Mursula et al., 2015). Accordingly, the correlation between geomagnetic activity and sunspots, especially in monthly averages as used here, is small, although significant (cc = 0.36,p < 10 −32 ).

Sudden Stratospheric Warmings
In this work we study how SSWs influence the dynamical response of the polar atmosphere to geomagnetic activity/EEP. SSW times are often identified (Charlton & Polvani, 2007;Holt et al., 2013) as times when the daily zonally averaged zonal wind at 60 • N and 10-hPa level turns easterly. In addition to this definition, for example, Charlton and Polvani (2007) (CP2007) also required that the zonal wind must have remained westerly at least 20 consecutive days before an SSW time and after the SSW recovery the zonal wind must remain westerly for at least 10 days. We slightly modified these criteria by requiring 10 consecutive days of westerly wind before SSW time and five consecutive days of westerly wind after the wind turns from easterly to westerly as the SSW recovers. This approach leads to only small differences in the number of identified SSWs compared to the CP2007 definition. Our identification criteria identify six additional late winter SSWs (compared to CP2007) where the zonal wind experiences large oscillation and momentarily becomes negative. Our criteria also identify one January SSW missed by the CP2007 definition because it has only 18 previous days of consecutive westerly wind instead of 20. These small differences compared to CP2007 do not significantly impact the results of this paper but do improve the statistical significance of the results slightly.
For each identified SSW, we also determined the month where the wind remains easterly longest and flagged the corresponding month as an SSW month. Typically, this is the month when the SSW starts, but sometimes, if the SSW begins near the end of the month, the SSW mainly affects the next month. Altogether, this process identified 42 SSWs.
Other previous studies (e.g., Maliniemi et al., 2013;Salminen et al., 2019) have found that when considering atmospheric responses to EEP or geomagnetic activity in monthly or winter-time averages using regression or correlation analysis, the obtained responses are influenced by the 2003/2004 winter, when an exceptionally strong SSW occurred on 5 January 2004. In these studies this particular winter is a clear outlier, which decreases the correlation between the polar vortex (or ground NAO index) and EEP. Also, in this study we find that this particular winter is a significant outlier and will be neglected from the first part of the analysis (section 5) that does not specifically consider the effects of SSWs. However, we will show that when carefully considering the effects of SSWs and their timing on the Ap-related responses in the atmosphere, the 2003/2004 SSW no longer constitutes an outlier.
We will also neglect in all analyses those winter months, which are potentially contaminated by effects of volcanic aerosols. These winter months were determined as those when the stratospheric aerosol optical depth averaged over the Northern Hemisphere (obtained from NASA/GISS at http://data.giss.nasa.gov/ modelforce/strataer/) exceeds 0.025 units. After excluding the winter months with volcanic contamination, we were left with 39 SSW events, out of which 4 occurred in December, 16 in January, 8 in February, and 11 in March. Note, however, that the results obtained in this paper are not greatly dependent on the removal of volcanic activity years and remain very similar even if all volcanic years are included or only the first year after volcanic eruption is excluded (as was done, e.g., by Salminen et al., 2019).

Methods
The effect of geomagnetic activity on the atmospheric parameters is studied here using linear regression analysis, which is performed separately in each latitude-pressure-level grid box. The regression is performed separately for four winter months (December, January, February, and March) in section 5 and for shorter time period averages associated with SSW times in section 6. The regression analysis method employed here has been used successfully in past studies and is described, for example, by Maliniemi et al. (2018). We review 10.1029/2019JD032137 here the main aspects of the methodology. Before regression, the atmospheric and Ap index data are first detrended by removing from the data a smooth trend estimated with the LOWESS method with a 31-year window (Cleveland & Devlin, 1988). This ensures that the results depict interannual variations rather than long-term trends. The regression is then performed using the Cochrane-Orcutt method (Cochrane & Orcutt, 1949), where the residual term is modeled as an autoregressive AR(1) process, instead of white noise as in regular regression. Thus, the regression model is where X t is the explaining variable (Ap index time series in our case), Y t is the time series of the considered atmospheric parameter in the given grid box, is constant, is the regression slope, and the residual term is t = t−1 + e t with being the lag-1 autocorrelation of the residual and e t being normally distributed white noise. The above regression model can also be written as The regression parameters ( and ) and can be solved iteratively by first performing normal regression using equation (11) and estimating from the residual time series t . Then, regression is performed using equation (12), and the value of is obtained in the previous step. New estimates for the parameters and can be obtained from the regression coefficients of equation (12). After that, a new estimate for can be obtained by calculating the new residuals t from equation (11). The procedure of estimating regression parameters from equation (12) and then can be iterated until converges. Because the residual term of equation (12) is uncorrelated white noise, the variances and statistical significance of the regression parameters of equation (12) can be estimated with a standard two-tailed Student's t-test (note that for equation 11, Student's t-test would be inappropriate due to the AR(1) noise). The Cochrane-Orcutt regression method is particularly suitable for problems where the model residual is not uncorrelated white noise but contains a physical signal from some other factors, which may have significant autocorrelation but are omitted as an explicit parameter from the model. Their effects can implicitly be modeled by the autoregressive residual term. Note that the AR(1) noise model is an idealization, which is often, but not always a good assumption. In order to test its validity, here we performed the significance analysis of the regression results also with a computationally much more time-consuming Monte-Carlo simulation. We generated random surrogate data for the Ap index and the atmospheric data by randomly mixing the phases of the Fourier components of each data series (Maliniemi et al., 2013;Thejll et al., 2003). This produces random time series with the same mean, variance, and autocorrelation functions as the original ones and furthermore maintains the spatial correlations in the atmospheric data. The regression fitting was then done for these surrogate data and repeated 10,000 times. The p values of the regression coefficients were then estimated with the histograms obtained from the Monte-Carlo repetitions as the fraction of values, which deviate from 0 more than the obtained regression coefficients. The p values obtained this way were very close to the p values from the Cochrane-Orcutt method, which validates the use of the AR(1) noise model here.
After obtaining the final regression coefficient values in each grid box, we then multiply them with the overall standard deviation of the monthly Ap index time series. The regression coefficients multiplied this way then yield the response in the atmospheric variable Y related to one standard deviation increase of monthly Ap index.
When considering solar-related effects on the stratosphere, a significant influence also comes from the solar UV irradiance, for which one often uses sunspot number as a proxy. In this study we do not discuss the effect of varying UV irradiance as it has been studied in many previous works (Gray et al., 2010). However, a concern may be raised about the validity of a regression analysis, which includes only the Ap time series. The validity of the results shown here was tested by including the sunspot number as an additional explanatory variable. Additional tests were also done by including El Niño/Southern Oscillation index, volcanic aerosols, and QBO as explanatory factors into the regression. The atmospheric responses to these additional factors will not be discussed in this paper, but in all of these tests, the Ap-related signals and their statistical significance remained rather similar to those reported here. This verifies the robustness of the results presented here. Similar results were recently obtained in a multilinear regression study of sea-level pressure and ground-level zonal winds (Maliniemi et al., 2018) indicating that the Ap-related effects are independent of solar UV effects and other internal climate drivers. This is also expected, since the Ap time series has only a rather weak correlation with SSN as mentioned above. Furthermore, the effect of the additional explanatory variables is implicitly modeled by the autoregressive residual employed here for the Ap-only regressions.
When performing the regression analysis for different winter months, we found that the signals are strongest in December and January if we do not use a lag for the Ap index. For February, the signals were found to be strongest if we use a lag of 1 month for the Ap index, and for March, the strongest signals were found with Ap lag of 3 months (although a lag of 2 months is fairly similar). These lags for Ap are thus employed in all analyses. Similar time lags were already previously found optimal (Maliniemi et al., 2013(Maliniemi et al., , 2014Salminen et al., 2019) and agree with the time it takes for majority of NO x /NO y created by energetic particle precipitation to descend from the mesosphere to the stratosphere (Funke et al., 2014a. Note that these lags also correspond well to the effective lags, which Funke et al. (2014b) found optimal for correlation between EPP produced NO y in the stratosphere/lower mesosphere and the Ap index. Figure 2 shows the response of monthly averaged zonal mean zonal wind and temperature in the Northern Hemisphere to one standard deviation increase in monthly averaged Ap in different winter months estimated with the regression method discussed in section 4. Responses similar to Figure 2 have already been reported in many previous observational (e.g., Salminen et al., 2019;Seppälä et al., 2013) and modeling studies (e.g., Arsenovic et al., 2016;Baumgaertner et al., 2011). In all winter months, one can see a significant Ap-related enhancement of westerly wind poleward of 40 • -50 • N and a corresponding weakening of the westerly wind at lower latitudes. These signatures extend from the stratosphere to the troposphere and even to ground level. The enhancement of westerly wind is associated with anomalous cooling of the lower polar stratosphere below about 50 hPa and anomalous warming of the upper polar stratosphere above 10 hPa. The temperature signals are not clearly visible in December, but starting from January, they become stronger and statistically significant (here, all winters are included, except the excluded volcanic years and the 2003/2004 SSW winter).

Response of Wintertime Polar Atmosphere to Geomagnetic Activity
The features seen in Figure 2 have been observed in earlier studies (e.g., Arsenovic et al., 2016;Baumgaertner et al., 2011;Christiansen et al., 1997), which suggest that the warming in the upper stratosphere may be caused largely by ozone loss-related reduction of radiative cooling. The lower altitude cooling has been suggested to result mostly from dynamical cooling (e.g., Baumgaertner et al., 2011) due to the weakening of downwelling within the polar vortex induced by changes in the propagation pattern of planetary waves (Limpasuvan et al., 2005). The temperature signals are in thermal wind shear balance with the westerly zonal wind so that the wind is enhanced between warming and cooling regions. The strengthened polar vortex also filters eastward propagating gravity waves, while allowing the propagation of waves with westward velocity into the upper stratosphere and mesosphere. When the waves break there, they impart easterly momentum on the flow, which leads to enhanced downwelling that further contributes to the heating in the upper Figure 3. The panels from top to bottom represent the response of zonally averaged zonal wind, temperature, integrated EP-flux divergence, and integrated adiabatic heating rate to one standard deviation increase of Ap index. Columns from left to right depict winter months from December to March, and the title of each column also indicates the number of data points used in the regression model for the corresponding month. The responses were calculated for those Decembers when an SSW occurs either in December or January (average date for those SSWs is 17 January) and for those January-March months when SSW occurred either in February or in March (average date for those SSWs is 6 March). Horizontal axes show the latitude in degrees, and vertical axes show the altitude in air-pressure units of hPa. The contours indicate the 90% (gray) and 95% (black) significance levels.
stratosphere by dynamic warming. Such changes in downwelling and dynamic heating have been reported, for example, by Limpasuvan et al. (2005).
We then proceed to study how the occurrence of SSWs during the winter affects the signals observed in Figure 2. Figure 3 displays results of the regression analyses for zonal wind, temperature, integrated EP-flux divergence, and integrated adiabatic heating rate for all SSW winters. For the analysis in December, we included winters when SSWs occur in either December or January and for January-March those winters when the SSWs occur in February or March.
The responses of zonal wind and temperature to the increased Ap index seen in Figure 3 are strong and significant in all winter months, even stronger than in Figure 2. They imply that geomagnetic activity is associated with enhanced polar vortex and associated bipolar temperature pattern (warming in the upper stratosphere and cooling in the lower stratosphere) especially in those winters when SSWs occur. The two bottom rows of Figure 3 display the corresponding responses in integrated EP-flux divergence and integrated adiabatic heating rate, respectively. One can see that in December-February, there is a strong and significant positive response in integrated EP-flux divergence in the middle stratosphere between 10 and 100 hPa poleward of 50-60 • N. There is also a significant negative signal at lower latitudes around 20 • -30 • N in the upper stratosphere between 1 and 10 hPa. Such statistically significant signals in EP-flux divergence are not seen in March. The bottom row of Figure 3 shows that there are also corresponding Ap-related signals in the integrated adiabatic heating rate. In December-February, there is significant adiabatic cooling in the lower polar stratosphere below about 30 hPa, which extends into troposphere in January and February. These cooling signals are located just below the signals of enhanced EP-flux divergence, which is consistent with the so-called downward control principle (Haynes et al., 1991), which states that downwelling at a certain altitude is determined by wave-driven zonal drag (EP-flux convergence or gravity wave drag) above that altitude. In March, one can also see statistically significant adiabatic warming in the upper polar stratosphere above 10 hPa. Note that the adiabatic heating rate in the polar stratosphere is mainly associated with changes in the downwelling of the polar air due to changes in the residual meridional circulation. The observed cooling (warming) signals are associated with corresponding weakening (strengthening) of the downwelling in the polar stratosphere (not shown explicitly).
Ozone loss is expected to cause radiative warming (reduction of long-wave cooling) in December-February and radiative cooling (reduction of short-wave warming) in March. It is interesting to note that the responses of zonal wind and temperature to the EEP in late winter, especially in March, are fairly similar to those in earlier winter months. This is evidently because the temperature response is greatly influenced by wave activity and resulting changes in meridional circulation. As seen in Figure 3 (bottom), the dominant warming of upper stratosphere in March is caused mainly by enhanced downwelling and subsequent dynamical heating so that it even dominates the expected ozone loss-related cooling in the upper stratosphere. The downwelling is likely driven by increased easterly gravity wave drag in the upper stratosphere/mesosphere, which is due to the enhanced polar vortex that filters westerly gravity waves. Note also that the March months in Figure 3 are taken from the same winters as the preceding Februaries so that the EEP-related enhancement of polar vortex in February still influences the wave-mean-flow interaction in subsequent March.
We also performed the analysis by selecting those winters, when an SSW occurs at any time during the winter (before or after each winter month under consideration), and the results were similar but clearly weaker than those in Figure 3. Note also that unlike in other recent studies (e.g., Maliniemi et al., 2013;Salminen et al., 2019), the January and February signals in Figure 3 are strong and significant even when including the exceptionally strong SSW that occurred in 5 January 2004 during a winter when the Ap index also had a very large value. These past studies found that this winter is a clear outlier, which tends to diminish the positive correlation between Ap index and polar vortex strength. In Figure 3, the January-March months of this particular winter do not contribute because the SSW must occur in February or March when considering January-March months. Figure 4 shows the corresponding atmospheric responses to Ap in those remaining winters/months that are not included in the analysis shown in Figure 3. In other words, Figure 4 contains those Decembers, when an SSW does not occur in December or January (but may or may not occur later in the winter) and those January-March months, when an SSW does not occur in February or March (but might have occurred before in the winter). One can clearly see that there are no strong or significant signals in polar region in these winters in December-February. In March, there is a statistically significant weakening of the stratospheric polar vortex, which is associated with warming of the polar stratosphere, enhanced planetary wave convergence, and enhanced downwelling leading to anomalous adiabatic warming.

Dynamical Amplification of the Ap Response Before SSWs
Overall, Figures 3 and 4 indicate that SSWs are necessary for the dynamical effect of Ap to be observed in the polar stratosphere. A possible explanation could be that the enhanced downwelling associated with mid-winter SSWs facilitates stronger descent of ozone-destroying NO x compounds created by EEP from the mesosphere into the stratosphere (e.g., Holt et al., 2013). This in turn would lead to a stronger EEP effect on the polar vortex sometime after the vortex has recovered from the SSW. However, the timing of the SSWs relative to the studied month in the above analysis suggests that the dynamical EEP-related effect takes place already before the SSW occurs (at least in December-February). If this is true, there must be a mechanism in the stratosphere, which not only leads to an SSW but also allows the EEP-related ozone loss to efficiently couple to the dynamics already before the SSW occurs. The responses observed in Figure 3 suggest that this mechanism may be related to planetary wave convergence and associated changes in residual circulation.
In order to study this in more detail, we considered the composite ERA-40/ERA-Interim data in daily resolution and computed averages of the different quantities in a 10-day time window from 15 to 5 days before each SSW time (including the January 2004 SSW). In addition, we calculated the average Ap index in the corresponding time window by applying the same time lags to the Ap as above in the monthly regression analysis (i.e., no lag for December or January, 1 month for February, and 3 months for March). Before calculating the averages, we removed the overall average seasonal variation from the atmospheric data. This was done by computing the overall average of the data in each grid box for each calendar month and interpolating the results for each day of the year (i.e., this seasonal variation is the same for all years). In addition,  Figure 3 but now using those remaining data points, which were not included in the regression analyses shown Figure 3. That is, the regressions were calculated for those Decembers where an SSW did not occur in December or January and those January-March months where an SSW did not occur in February or March.
we removed from the Ap index and atmospheric data the 31-year smooth trend (see section 4). Figure 5 shows the results of the regression analysis for these 10-day averages of zonal wind, temperature, integrated EP-flux divergence, and integrated adiabatic heating rate. The responses are quite similar to those for December-February in Figure 3. They indicate that enhanced geomagnetic activity is associated with a strong polar vortex enhancement before the SSW. The polar vortex enhancement is accompanied by significant anomalous planetary wave divergence in the stratosphere somewhat below the region of strongest zonal wind enhancement. The anomalous planetary wave divergence causes relative acceleration of westerly zonal wind according to equation (6) and weakening of downwelling, which is indicated by the significant adiabatic cooling anomaly in the lower stratosphere. The adiabatic cooling leads to the negative temperature response in the lower polar stratosphere and further enhances the westerly zonal wind according to the thermal wind shear balance. When studying these responses related to the SSW times, we also tried to extend the time window backward and forward in time. We found that especially, when extending the window forward in time closer to the SSW time (or after the SSW), the Ap-related signals quickly diminished, indicating that the Ap-related response is best observed slightly before the SSW.
The above analyses indicate that the effect of geomagnetic activity (or EEP) on the polar vortex is related to the strong wave-mean-flow interaction, which amplifies the initial vortex enhancement related to EEP-induced ozone loss to the observed level by planetary wave refraction. This dynamical effect seems to be observable slightly before the SSW occurs. In order to understand why the pre-SSW times are important for the mechanism to take place, let us next consider how the atmospheric conditions during these periods differ from other times. Figure 6 shows the average differences between pre-SSW times (5-15 days before SSW central times) and those full winters (December-March) when there are no SSWs. Statistical significance of the results was estimated by Monte-Carlo bootstrapping with 1,000 repetitions. For each repetition, we first randomly generated in daily resolution the same number of artificial SSW times as we have true SSW times. For each such time, we took the data from 5 to 15 days before and finally computed their average. We also computed the average of the atmospheric parameters in all those full winters (December-March), which did not contain the randomly generated artificial SSW times. The statistical significance of the observed difference between pre-SSW times and other winter times was then estimated by determining the fraction of the 1,000 repetitions where the random differences are larger than the observed ones. From Figure 6, one can see that the pre-SSW times are characterized by weaker than average stratospheric zonal wind in the latitudes above 40 • N. Figure 6 also shows that pre-SSW times are characterized by strong planetary wave convergence throughout the stratosphere above 100 hPa poleward of 40-50 • N latitude. This wave convergence drives enhanced downwelling in the polar region, which is seen as enhanced adiabatic heating and increased polar stratospheric temperature. The enhanced downwelling/warming in the polar region is accompanied by a weaker signal of enhanced upwelling/cooling in the lower latitudes equatorward of 30 • N.
The overall state of the stratosphere during pre-SSW times is associated with enhanced planetary wave activity from the troposphere, which is enhanced by weakened zonal wind that allows further planetary wave activity to converge on the polar vortex (Matsuno, 1970). During non-SSW times, the wave activity into the stratosphere is much weaker. It thus appears that the reason why the Ap/EEP effect on the polar vortex is observed during pre-SSW times is the strong planetary wave activity, which can dynamically amplify the enhancement of the polar vortex due to radiative changes related to EEP-induced ozone loss. Conversely, it seems that if a sufficient amount of planetary wave activity is not present, the initial polar vortex enhancement related to EEP-induced ozone loss is not strong enough to be distinguished from internal variability of the polar vortex. It should be noted that amplification of the response is observed in the superposed epoch analysis above most clearly during the 5 to 15-day period before the SSWs because wave activity is then systematically elevated in each SSW. Earlier than this period, wave activity varies greatly from one SSW event to the next, and no systematic signal is observed. In the monthly averages shown in Figure 3, the responses are better seen because each month includes some periods of time, which are closely preceding some SSWs. Figure 6. Differences of zonal wind, temperature, integrated EP-flux divergence, and integrated adiabatic heating rate between pre-SSW times (5 to 15 days before SSW) and those full winters when an SSW does not occur. Horizontal axes show the latitude in degrees, and vertical axes show the altitude in air-pressure units of hPa. The contours indicate the 90% (gray) and 95% (black) significance levels estimated by Monte-Carlo bootstrapping with 1,000 repetitions.
The dynamical amplification of the EEP response is a result of the wave-mean-flow interaction. It is known that planetary waves tend to be refracted away from stronger westerly winds and toward weaker westerly winds (e.g., Matsuno, 1970). Thus, if EEP initially, by ozone loss, influences radiative heating rates so that the westerly zonal wind is slightly enhanced by thermal wind shear balance, this change tends to enhance the divergence of planetary waves away from the region of zonal wind enhancement. This enhancement of EP-flux divergence will then accelerate the westerly wind further (see equation 6) and lead to further refraction of planetary waves, forming a positive feedback loop. It is also interesting to note that since the wave divergence is associated with weakened downwelling, the transport of ozone into the lower stratosphere is also weakened. This should result in EEP-induced dynamical reduction of ozone in the lower polar stratosphere, which would enhance the net ozone reduction on top of the initial chemical ozone loss due to NO x /NO y . Such dynamical ozone reduction was also suggested recently by Salminen et al. (2019). The initial EEP-related chemical ozone loss is not expected to occur in the middle or lower stratosphere especially in December-January but is most likely confined to the mesosphere and upper stratosphere (above 10 hPa), where the downwelling NO x /NO y reaches already in December . Thus, it is the initial EEP-induced changes in the mesosphere/upper stratosphere that initiate the suggested wave amplification mechanism in the presence of suitably strong wave activity. The resulting dynamical changes relatively quickly propagate down to lower stratosphere and even troposphere by means of wave-mean-flow interaction. In monthly averages or even 10-day averages, these changes would seem to be instantaneous, since the thermal relaxation times at these altitudes are less than 1 month (Mlynczak et al., 1999) and propagation time of planetary waves between ground and stratosphere is less than a week (e.g., Perlwitz & Graf, 2001).
In addition to enhancing the polar vortex when EEP is high, the wave-mean-flow interaction can also be thought to work in the opposite direction when the EEP is very low. In this case EEP-associated ozone loss is smaller than average, and the polar vortex initially tends to be slightly weaker than average. The weaker

Discussion and Conclusions
In this paper we have studied the response of the polar stratosphere to geomagnetic Ap index, which was used as a proxy of energetic particle precipitation. In agreement with previous studies we confirmed that an increase of EEP is associated with polar vortex enhancement and a corresponding warming (cooling) of the upper (lower) polar stratosphere in all winter months from December to March. We then showed that the EEP-related response in the atmosphere arises from those winters when an SSW occurs. Conversely, in winters without an SSW, zonal winds do not display any appreciable response to EEP. Interestingly, the EEP-related signal was found to take place before the SSW occurs, thereby indicating that certain atmospheric conditions preceding the SSW allow the EEP effect to be observed.
We also studied the corresponding EEP responses in planetary waves and adiabatic heating anomalies, which are also related to changes of meridional circulation. We found that at times slightly preceding the SSWs, the EEP response is characterized by clear and statistically significant enhancement of planetary wave divergence in the stratosphere. This increase of divergence leads to a relative acceleration of the westerly zonal wind and to a weakening of downwelling in the polar stratosphere, which causes adiabatic cooling of the lower polar stratosphere. Cooling further enhances the westerly zonal wind by thermal wind shear balance. These observations clearly indicate that planetary wave refraction and corresponding changes in downwelling in the polar stratosphere are key ingredients in generating the observed polar vortex enhancement as a response to increased EEP. It is likely that, by this wave-mean-flow interaction, the small initial changes in heating rates related to EEP-induced ozone loss are dynamically amplified to a level high enough to be distinguished from internal variability.
We then considered the atmospheric conditions preceding the SSWs that allow the above-described dynamical amplification of the EEP response to take place. We found that the pre-SSW times are characterized by a weaker than average polar vortex and strong convergence of planetary wave activity in the middle and upper stratosphere at mid-to high latitudes. Strong planetary wave convergence also drives enhanced meridional circulation throughout the stratosphere, which is seen as stronger downwelling (upwelling) and consequent adiabatic warming (cooling) in the polar (subtropical) stratosphere. It is likely that the stronger planetary wave convergence in the stratosphere, which is most systematically observed during the pre-SSW periods, allows stronger wave-mean-flow interaction to take place and thereby allows the EEP responses to be dynamically amplified. During times that do not lead to SSWs, the planetary wave activity is apparently not strong enough to amplify the EEP response above internal variability.
Several past studies have shown that mid-winter SSWs enhance the downwelling of EEP-created NO x /NO y from the mesosphere into the stratosphere and thus enhance the EEP-related ozone loss sometime after the polar vortex has recovered (e.g., Funke et al., 2016;Holt et al., 2013;Päivärinta et al., 2016). Our results show that for the polar vortex dynamics, the pre-SSW times are perhaps even more important than the enhanced period of downwelling when the SSW has already happened. In order to better understand the impact of energetic particle precipitation on atmospheric dynamics and not only on atmospheric chemistry, one needs to consider the pre-SSW times, which are characterized by strong wave-mean-flow interactions driven by planetary waves. This is also true for modeling studies and may be important to take into account if the model is not capable of producing realistic SSWs at correct times relative to the energetic particle input.
There is another interesting implication of our results. If sufficiently strong planetary wave activity is needed to amplify the EEP response, then any factor internal or external to the climate system, which significantly influences planetary wave activity into the stratosphere, can possibly act as a modulator for the energetic particle effect. One such well-known factor is the QBO, which in its easterly phase is known to guide more planetary wave activity into the polar stratosphere (Holton & Tan, 1980). Many past studies have reported that the EEP-related effects on the winter time polar vortex and associated ground climate variations are stronger in the easterly QBO phase (e.g., Maliniemi et al., 2013Maliniemi et al., , 2016Palamara & Bryant, 2004). Salminen et al. (2019) suggested that a possible explanation for the QBO modulation is the stronger meridional circulation and associated planetary wave activity, which makes conditions more favorable for stronger wave-mean-flow interaction. Our results here give strong evidence in favor of this hypothesis. Internal stratospheric conditions can also strongly modulate planetary wave refraction even if wave activity impinging