Evidence for stratospheric sudden warming effects on the upper thermosphere derived from satellite orbital decay data during 1967–2013

We investigate possible impact of stratospheric sudden warmings (SSWs) on the thermosphere by using long‐term data of the global average thermospheric total mass density derived from satellite orbital drag during 1967–2013. Residuals are analyzed between the data and empirical Global Average Mass Density Model (GAMDM) that takes into account density variability due to solar activity, season, geomagnetic activity, and long‐term trend. A superposed epoch analysis of 37 SSW events reveals a density reduction of 3–7% at 250–575 km around the time of maximum polar vortex weakening. The relative density perturbation is found to be greater at higher altitudes. The temperature perturbation is estimated to be −7.0 K at 400 km. We show that the density reduction can arise from enhanced wave forcing from the lower atmosphere.


Introduction
A stratospheric sudden warming (SSW) is a large-scale meteorological disturbance that usually takes place in the arctic region during winter [e.g., Andrews et al., 1987]. Upward propagation of planetary waves transfers energy and momentum from the troposphere to the stratosphere to initiate the disturbance [Matsuno, 1971]. An SSW significantly alters the mean state of the middle atmosphere [Liu and Roble, 2002], which in turn affects the propagation of other atmospheric waves into the upper atmosphere [Stening et al., 1997;Liu et al., 2010].
Ionospheric effects during SSW events have been recognized in recent years, especially in the wake of an unusually strong and prolonged SSW event in January 2009 [Manney et al., 2009]. Studies found marked changes in the plasma density at low latitudes [e.g., Chau et al., 2010;Goncharenko et al., 2010]. These ionospheric modulations have been attributed mainly to electrodynamic effects induced by enhanced semidiurnal tidal waves [e.g., Jin et al., 2012;Wang et al., 2014].
In contrast to the successful progress in the understanding of ionospheric effects during SSWs, it has been controversial whether an SSW has any measurable impact on the thermosphere. Unlike the ionosphere, remote sensing of the thermosphere is difficult, and available data are very limited. Using satellite-borne accelerometer data, Liu et al. [2011] suggested that the thermospheric density at 325-475 km underwent an anomalous decrease of 30-45% during the January 2009 SSW. This claim was, however, immediately questioned by Fuller-Rowell et al. [2011], who showed that the apparent density reduction reported by Liu et al. [2011] can be largely explained by changes in geomagnetic activity. It was therefore concluded that there was no evidence of thermospheric effects during the January 2009 SSW. Conde and Nicolls [2010], examining temperature data at 240 km for the same event, also noted difficulties in separating contributions of the SSW from other sources.
Another issue regarding possible SSW impact on the thermosphere is the lack of a consensus among predictions by general circulation models. Pedatella et al. [2014a] compared simulation results for the January 2009 SSW event obtained from these whole atmosphere models: GAIA [Jin et al., 2011], HAMMONIA [Schmidt et al., 2006], WAM [Fuller-Rowell et al., 2010], and WACCM-X . The comparison revealed that although the models show reasonable agreement in the troposphere and stratosphere, the model results become diverse as the altitude increases and have significant discrepancies in the thermosphere. Therefore, there is no reliable prediction at present about how the thermosphere would respond to SSWs and how much the effect would be.
In this study, we use long-term records of the global average thermospheric density derived from satellite orbital decay data, which for the first time make it possible to investigate the thermospheric response to SSW events in a statistical way. The results reveal a global reduction in upper thermospheric density during SSWs. We suggest enhanced wave forcing from below as a possible cause for the density reduction. We use the thermosphere-ionosphere-electrodynamics general circulation model (TIE-GCM) [Richmond et al., 1992] to demonstrate how wave forcing would affect the upper thermosphere.

Data and Models
The height-dependent, global average density data set we use was detailed by Emmert [2009Emmert [ , 2015. Briefly, the decay rates of the mean orbital period of near-Earth orbit satellites are proportional to the total mass density along the orbital track. This relationship can be used to evaluate the thermospheric density from orbital trajectory data. The use of historical orbital records enables construction of long time series. Combining densities from multiple objects, it is possible to obtain a representation of the global average density, but with a limited temporal resolution. The data set constructed by Emmert [2015] incorporates orbital data from ∼5000 objects. It covers the period from January 1967 to December 2013 and altitudes between 250 and 575 km. Estimated daily relative accuracy is ∼2%, and estimated absolute accuracy is ∼10% [Emmert, 2009, section 4]. Although we use daily data, the maximum temporal resolution is 3 days because of 3 day cubic spline smoothing involved in the density retrieval procedure [Emmert, 2009, section 7]. The reader is referred to Emmert [2009] for further information on the derivation of the global average mass density and error evaluation.
The Global Average Mass Density Model (GAMDM) is an empirical model of global daily average thermospheric density, which was constructed on the basis of function fits to the density product described in the previous paragraph. Emmert and Picone [2010] discussed the construction details of GAMDM. The model uses solar activity index F 10.7 , geomagnetic activity index Kp, and day of year as input parameters to evaluate contributions of solar and geomagnetic activity, and season. The residuals between the data and GAMDM may be interpreted as resulting from other sources, such as increasing CO 2 concentration [Emmert et al., 2008]. The latest version of GAMDM, which we use, takes into account the effect of the CO 2 increase [Emmert et al., 2014;Emmert, 2015]. Since GAMDM does not include any forcing term associated with SSWs, the effect of SSWs on the thermospheric density, if there is any, would lead to overestimation or underestimation of the model prediction in comparison with the data. As reported by Emmert et al. [2014], the density data during 2006-2010 were consistently lower than GAMDM predictions, by up to 20% on average at 400 km [Emmert, 2015, Figure 2c]. In addition, the 400 km residuals have an interannual variance of ∼3% [Emmert et al., 2014, Table 5]. Since these annual-scale anomalies are presumably unrelated to shorter-term SSW events, we remove them by subtracting running annual averages of the residuals from the residual time series prior to conducting our analysis.
Different studies have used different criteria for the detection of SSWs [Butler et al., 2015]. One of the most commonly used diagnostics is the zonal mean zonal wind at 60 ∘ and 10 hPa (∼32 km), which has been found useful in studies of stratosphere-troposphere coupling [e.g., Charlton and Polvani, 2007]. For upper atmospheric studies, Zhang and Forbes [2014] and Chau et al. [2015] used the zonal mean zonal wind at 70 ∘ at 1 hPa (∼48 km) to find a remarkable correlation between the polar vortex weakening and semidiurnal lunar tidal activity in the lower thermosphere, in terms of both timing and intensity. The success of the use of zonal wind at 1 hPa, instead of 10 hPa, probably reflects the fact that a breakdown or weakening of the polar vortex begins at an altitude higher than 10 hPa and gradually descends to lower levels. In this study, we also use the zonal mean zonal wind at 70 ∘ at 1 hPa to identify SSWs. Wind and temperature data were obtained from global meteorological reanalysis produced by the European Centre for Medium-Range Weather Forecasts (ECMWF). Specifically, we use ERA-Interim [Dee et al., 2011] for 1979and ERA-40 [Uppala et al., 2005 for 1966-1978.
The TIE-GCM is a first-principles model of the coupled thermosphere-ionosphere system [Richmond et al., 1992;Qian et al., 2014]. It solves the continuity, momentum, and energy equations with self-consistent electrodynamics over a height range from approximately 96 km to 600 km. For the present study, we use version 1.95 of the TIE-GCM. The model resolution is 5 ∘ in latitude and 5 ∘ in longitude with two grid points per scale height. Tidal forcing at the lower boundary was specified for migrating diurnal and semidiurnal tides using the Global Scale Wave Model (GSWM) Forbes, 2002, 2003]. Also, eddy diffusivity at the lower boundary was specified in the way described by Qian et al. [2009].

Case Study: January 2009 SSW Event
We first present the thermospheric density variability during the January 2009 SSW event. Figure 1a depicts the zonal mean zonal wind at 70 ∘ N at 1 hPa. The solid line shows the wind for the winter period of 2008-2009, while the dashed line shows the climatological median derived from the daily data during 1966-2014. A reduction of the zonal mean zonal wind from the climatological median is evident around day 20-30, indicating significant weakening of the polar vortex, and hence the occurrence of SSW. We plot in Figure 1b the thermospheric density at 400 km for the corresponding period. The results from the data and GAMDM are shown by the black and blue lines, respectively. It can be seen that GAMDM captures most of the density variability. Some of the changes are clearly correlated with changes in geomagnetic activity, which is shown in Figure 1c. GAMDM overestimates the density around the SSW event, especially after the peak polar vortex weakening. The maximum discrepancy occurred on day 25, when the data-GAMDM residual was −7% (after the annual average was removed from the residual time series).
It is obviously arguable if the density reduction around day 25 is actually associated with the SSW event. One can see in Figure 1b that GAMDM similarly overestimates the density around day −40 and day −30 without SSW. According to Emmert et al. [2014, Table 5], the day-to-day variability (standard deviation) of the GAMDM residuals is ∼13% at 400 km. Therefore, the density reduction of 7% is smaller than the background noise.

10.1002/2015GL065395
Our investigation of other individual SSW events also encountered difficulties in separating possible SSW signals from the large noise. We thus choose not to conclude about SSW effects on the thermospheric density from individual events. Instead, we take a statistical approach that reduces the noise, as described in the next section.

Superposed Epoch Analysis
A superposed epoch analysis is a statistical analysis technique, widely used in geophysical research, to estimate the characteristic response of a system to a type of event [e.g., Fejer et al., 2002;Bristow and Jensen, 2007]. We employed the technique to estimate the average response of thermospheric density perturbation to SSW events during 1967-2013. An SSW event was identified when the daily mean value of the zonal mean zonal wind at 70 ∘ and 1 hPa was below the climatological median by more than 20 m/s for at least three consecutive days. We focused on the Northern Hemisphere mid winter period, defined here as ±45 days from the first day of the year. Once an SSW event was detected, the day for the peak polar vortex weakening (i.e., minimum zonal mean zonal wind) was assigned as "day 0," and no event was considered for the following 30 days. We found 37 SSW events during the period examined, which leads to the occurrence rate of 0.79 events per winter, in agreement with the results from other SSW definitions in the literature [e.g., Butler et al., 2015]. A list of the 37 SSW events is given in the supporting information. Figure 2a shows the average of the wind anomaly for the 37 SSW events. The westward wind anomaly starts to develop 5-10 days before the peak polar vortex weakening. As the altitude increases, the magnitude of the westward wind anomaly becomes greater, and the peak of polar vortex weakening occurs earlier. After the peak polar vortex weakening, the wind gradually recovers to the normal level in 20-25 days at 1 hPa.
In Figures 2b-2e, the black lines depict the average of the density anomaly for the 37 SSW events at 275, 375, 475, and 575 km. Similar plots showing the density anomaly for all 37 individual SSW events can be found in the supporting information. The average signal reveals a reduction of the thermospheric density during SSWs at all heights. The maximum density reduction occurs 1-4 days after the peak polar vortex weakening. The relative density perturbation during SSWs is greater at higher altitudes.
Figures 2f-2i display the distribution of the density anomaly for days 1-4. The average density perturbations are −2.9%, −4.4%, −5.9%, and −6.9% at 275, 375, 475, and 575 km, respectively. The density anomaly has an approximately Gaussian distribution, especially at 475 and 575 km. The 1 uncertainty of the mean may be estimated to be the standard error of the mean (SEM), i.e., the standard deviation of the sample divided by the square root of the number of the independent samples (i.e., 37). The calculated 1 uncertainty estimates are 1.0%, 1.5%, 2.0%, and 2.4% at 275, 375, 475, and 575 km, respectively. It is noted that these uncertainty estimates relate to the mean for the 37 SSW events, not individual events. In order to ensure that these uncertainty estimates are reasonable, the 1 uncertainty of the mean was also evaluated using a Monte Carlo simulation. We selected 37 random epochs in the data set and calculated the average density anomaly. This procedure was repeated 2000 times. The standard deviation of these 2000 mean values was computed as the estimated 1 uncertainty. It is noted that the standard deviation of the data-GAMDM residuals do not depend on solar activity, season, or geomagnetic activity [Emmert and Picone, 2010, Figure 10]. The derived 1 uncertainties are 1.3%, 1.8%, 2.4%, and 2.9%, which is largely consistent with the SEMs. Thus, the deviation of the average density perturbation during SSWs is larger than the estimated 2 uncertainty of the mean at all heights.
We also compared the average density perturbation for the first 19 SSW events with that for the last 18 events, which are plotted in Figures 2b-2e by the green and orange lines, respectively. The pattern is remarkably consistent between days −15 and 7, with the maximum negative response near the peak polar vortex weakening. The consistency between the two halves of the data set gives us further confidence that the density reduction during SSWs is physically meaningful.

Temperature Response
Simulations and observations have indicated significant temperature changes in the polar middle atmosphere during SSWs [e.g., Liu and Roble, 2002;Funke et al., 2010]. These studies showed mesospheric cooling (70-90 km) and lower thermospheric warming (120-140 km), along with stratospheric warming (30-50 km). However, little is known about how the thermospheric temperature changes on a global scale during SSWs. Liu et al. [2013 reported a long-term decrease (∼30 days) in the global average temperature by −12 K for the January 2009 event using GAIA. There has been no observational evidence to verify the results. . Height dependence of the thermospheric total mass density anomaly during SSWs averaged over days 1-4 (blue crosses). Results for the TIE-GCM numerical experiment are also shown; changes in the global mean total mass density (magenta dashed line), temperature (green solid line), and one half the plasma density (black dash-dotted line) due to elevated eddy diffusivity at the lower boundary. The plasma density values were scaled for better display.
Here, we estimate the change in the global average temperature during SSWs from the height gradient of the density perturbation. Under hydrostatic equilibrium, the mass density is related to temperature T and mean molecular mass M by where z is height, g is the acceleration due to gravity, and k is Boltzmann constant (= 1.38 ×10 −23 J/K). At 400 km, the gradients in T and M are small and can be neglected. Also ignoring perturbations in M, the perturbed equation is In Figure  Thus, ΔT is −7.0 ± 2.5 K, where the 1 uncertainty was derived from the estimated 1 uncertainty in the vertical gradient of the density perturbation. The actual 1 error in the temperature perturbation may be greater as we did not take into account errors in T and M. Our estimate for the temperature perturbation is similar to but somewhat smaller than the 12 K reduction predicted by Liu et al. [2013 for the January 2009 event using GAIA.

Mechanism for Density Reduction
Changes in the thermospheric density at a fixed altitude can arise from changes in both temperature and composition [e.g., Lei et al., 2010]. The mechanism we focus on in this study is wave forcing from below the thermosphere. The dissipation of upward propagating gravity waves in the mesosphere and lower thermosphere induces a downward transfer of heat, which causes cooling at the higher levels [Walterscheid, 1981;Akmaev, 2007]. Cooling at a given height of the thermosphere would lead to a reduction of the density in the region above owing to the contraction of the air. Besides, breakdown of gravity waves generates diffusive turbulence [Lindzen, 1981], which induces a downward transport of atomic oxygen O. Since the atomic oxygen is the major constituent of the thermosphere at 250-575 km, the loss of O due to the downward transport means a reduction of the total mass density. Figure 3 illustrates the effect of enhanced gravity wave dissipation in the lower thermosphere. We performed two TIE-GCM simulations. A "Base Case" simulation was run for steady state 1 January conditions with F 10.7 = 135 sfu and Ap= 12 nT, and the other simulation was run with the same model configuration except that the eddy diffusivity at the model lower boundary (∼96 km) was increased by 20%. We calculated the global mean of total mass density, temperature, and plasma density for each run. Plotted in Figure 3 are the differences. The modulation of the eddy diffusivity at the TIE-GCM lower boundary has been used in previous studies to evaluate the effect of gravity wave forcing on the thermosphere-ionosphere system [e.g., Qian et al., 2009Qian et al., , 2013. The results in Figure 3 reveal a reduction in the global mean density, similar to observations during SSWs. It is known that changes in the background wind during SSWs significantly influence upward propagation of gravity waves. Studies have shown enhanced gravity wave forcing during the developing phase of SSWs [e.g., Hoffmann et al., 2007;Yamashita et al., 2010;Yigit and Medvedev, 2012].
Although it is not known how the eddy diffusivity changes in the mesosphere and lower thermosphere during SSWs, the short-term change of 20% seems plausible given that the eddy diffusivity varies (locally, at least) by a factor of 3 or so on the seasonal time scale [Kirchhoff and Clemesha, 1983;Fukao et al., 1994;Sasi and Vijayan, 2001].
SSWs also affect upward propagation of various tides and planetary waves [e.g., Pedatella et al., 2014a]. Particularly, the solar migrating semidiurnal tide has been reported to exhibit significant amplification during SSWs [e.g., Wang et al., 2012;Pedatella et al., 2014b]. The dissipation of upward propagating tides in the mesosphere and lower thermosphere is known to have an effect similar to that achieved by increasing eddy diffusivity [e.g., Akmaev and Shved, 1980;Forbes et al., 1993]. Our TIE-GCM experiment revealed that the results similar to Figure 3 can be obtained by increasing the amplitude of the solar migrating semidiurnal tide at the lower boundary by a factor of 3 (not shown). TIE-GCM experiments of this kind have been used in previous studies for the assessment of the thermospheric response to tidal forcing from below [e.g., Yamazaki and Richmond, 2013;Jones et al., 2014].
As shown in Figures 3, the enhanced wave forcing brings about a decrease in the F region plasma density. This is because the loss of O reduces the production of O + that dominates the F region plasma population. Satellite observations during the January 2009 SSW revealed a global-scale reduction in the electron density [Pancheva and Mukhtarov, 2011;Lin et al., 2012].
Although we have shown that enhanced wave forcing can induce a reduction of the global mean thermospheric density similar to observations during SSWs, it is open to question what waves would play a significant role and how important this mechanism would be relative to other mechanisms. Other possible mechanisms include cooling and shrinking of the mesosphere, which would lead to a reduction of the thermospheric density at a fixed altitude. At high latitudes, SSWs are often accompanied by cooling in the mesosphere. However, the region of the strong cooling is mostly confined to the Arctic region and the impact on the global mean density is expected to be marginal. Cooling/heating in the tropical mesosphere during SSWs is not well established, and its possible impact on the global mean thermospheric density is yet to be studied.