Role of gravity waves in vertical coupling during sudden stratospheric warmings

Gravity waves are primarily generated in the lower atmosphere, and can reach thermospheric heights in the course of their propagation. This paper reviews the recent progress in understanding the role of gravity waves in vertical coupling during sudden stratospheric warmings. Modeling of gravity wave effects is briefly reviewed, and the recent developments in the field are presented. Then, the impact of these waves on the general circulation of the upper atmosphere is outlined. Finally, the role of gravity waves in vertical coupling between the lower and the upper atmosphere is discussed in the context of sudden stratospheric warmings.


INTRODUCTION
The lower atmosphere, where meteorological processes take place, is the primary source of internal atmospheric waves: gravity waves (GWs), planetary (Rossby) waves, and solar tides. These waves can propagate upward and influence the dynamics and thermal state of the middle and upper atmosphere (see e.g., the reviews of Fritts & Alexander 2003;Laštovička 2006;Yigit & Medvedev 2015). Waves transfer their energy and momentum to the mean flow via breaking and dissipative processes, such as radiative damping, eddy viscosity, nonlinear diffusion, molecular diffusion and thermal conduction, and ion drag (Yigit et al. 2008). Sudden stratospheric warmings (SSWs) are spectacular events that disturb the circulation in the winter hemisphere. They affect not only the stratosphere, but their influence extends to the mesosphere and thermosphere. In the upper atmosphere, plasma processes, such as Joule and auroral heating, ion friction are important processes that shape the morphology and dynamics. Thus, interactions between the lower and upper atmosphere should be considered within the framework of the atmosphere-ionosphere system.
Such coupled upper atmosphere-ionosphere system is subject to the following internal and external influences: • Meteorological effects that encompass internal wave impacts and transient processes of lower atmospheric origin, • Internal processes due to nonlinearity, • Space weather effects that are associated with the solar and magnetospheric phenomena.
Among the meteorological effects, we distinguish a direct influence of internal GWs on the upper regions of the atmosphere. Although transient events such as SSWs are technically categorized as stratospheric processes, and, thus take place above the region of weather-dominated phenomena, they are often referred to as meteorological effects in the context of the upper atmosphere research.
The thermosphere-ionosphere system is highly nonlinear. In the real atmosphere, ion and neutral parameters vary simultaneously, and the resulting changes in the heating ought to contain higher order terms, which is indicative of the nonlinear nature of the system (Yigit & Ridley 2011a). The atmosphere-ionosphere system is subject to the influence of space weather, which can enhance these nonlinear processes and impact the upper atmosphere (Prölss 2011, and references therein).
In this paper, we report on the recent advances in understanding the meteorological effects in the upper atmosphere, focusing primarily on the links between SSWs, small-scale GWs, and thermosphere-ionosphere dynamics.

INTERNAL GRAVITY WAVES
Internal gravity waves are characteristic features of all stably stratified planetary atmospheres. GWs in the upper atmosphere have been studied for more than 50 years since the early work of Hines (1960). Their importance for the general circulation of the middle atmosphere has been greatly appreciated (e.g., Garcia & Solomon 1985;Becker 2011). However, despite the previous theoretical approaches to GW propagation into the thermosphere (Klostermeyer 1972;Hickey & Cole 1988), only since recently, the role of GWs in coupling the lower and upper atmosphere is being increasingly acknowledged Vadas & Liu 2009;Fritts & Lund 2011;Yigit et al. 2012a;Miyoshi et al. 2014;Hickey et al. 2011Hickey et al. , 2010Heale et al. 2014).
GWs are always present in the lower and upper atmosphere, however their amplitudes and dynamical importance differ with height. Wave energy is proportional to air density, and, therefore, a conservatively propagating harmonic has a larger amplitude in regions with lower density. In the troposphere, GW amplitudes are relatively small, however, their dynamical importance increases with height, and can no longer be neglected in the middle and upper atmosphere.
We next discuss basic principles of how GW processes are represented in atmospheric models, reviewing the underlying assumptions and limitations.
2.1. Principles of parameterization of gravity wave processes in global atmosphere models Spatial scales of GWs are considerably smaller than the planetary radius. Their sources are highly intermittent, and propagation is strongly dispersive. Therefore, the GW field in the thermosphere is highly irregular and transient. Unlike with distinct larger-scale planetary waves, it appears as an ever changing "sea of waves" with occasional well-defined and detectable packets. In many applications, such chaotic wave field and its influence on the larger-scale flow can be conveniently described in terms of statistical quantities devoid of the phase information. Examples of the most widely used statistical characteristics for the GW field are the variance φ 2 , vertical flux of horizontal momentum u w , sensible heat flux w T , etc, where w , T , and φ are the deviations of vertical velocity, temperature and of any field variable from the corresponding mean values, respectively.
General circulation models (GCMs) have spatial resolutions usually much coarser than the scales of GWs. Only few GCMs have endeavored to perform simulations with grids small enough in an attempt to resolve at least a part of the GW spectrum (e.g., Miyoshi & Fujiwara 2008;Miyoshi et al. 2014). In most simulation studies, the effects of subgrid-scale GWs have to be parameterized. This means that (1) the average effects must be presented in terms of statistical quantities similar to the described above, and the quantities have to be functions of the background flow. In other words, the parameterization has to self-consistently capture responses of the wave field to the evolution of the resolved larger-scale flow.
(2) Parameterizations should preferably be based on first principles, that is, they should rely on rigorous laws of physics rather than on a set of empirically introduced (tuning) parameters. Obviously, no parameterization can be devoid of such parameters as they are a substitute for an unknown. But the lesser the number of tunable parameters, the more sophisticated the parameterization is.
(3) Parameterizations must be verifiable. This condition means that they have to provide quantities, which can be compared with observations. For instance, GW-induced heating/cooling rates are hard to measure, but temperature variances T 2 can be.

Assumptions and limitations in gravity wave
parameterizations In modeling, it is assumed that the majority of GWs are generated in the lower atmosphere. Amplitudes of those excited in the upper layers and propagating downward decrease exponentially with height together with their influence on the mean flow. Therefore, (1) only harmonics propagating upward are considered in parameterizations. This assumption allows one to omit a detailed consideration of the wave reflection, and to (2) apply the Wentzel-Kramers-Brillouin (WKB) approximation. Under the WKB method, (3) only those harmonics are considered whose vertical wavelengths are much shorter than vertical variations of the background fields. Mathematically, the latter can be expressed as k z H 1, where k z is the vertical wavenumber and H is the density scale height. This limitation becomes very restrictive in the thermosphere, because fast (and long vertical-wavelength) harmonics have more chances to penetrate from tropospheric heights. In the real world, GWs propagate obliquely with respect to the surface. However, because k z k h for most harmonics, k h being the horizontal wavenumber, parameterizations (4) usually assume vertical-only propagation. Limitations of this approximation in the middle atmosphere have been recently discussed in the work by Kalisch et al. (2014), and higher-order effects have been found with a scheme employing ray tracing (Song et al. 2007). A special care should be taken with parameterizations extending to the thermosphere, where longer vertical wavelength harmonics (lower k z ) tend to propagate to from below. In other words, all gravity waves accounted for by a parameterization must remain within their grid columns. Finally, (5) all column-based parameterizations employ a steady state approximation. That is, transient processes of wave propagation assume an instantaneous response to changes in the forcing below. This approximation is suitable for modeling the general circulation, however, implications of time delay due to the finite group speed of wave packets should be carefully weighted for simulations of more rapid processes.
Parameterizations compute vertical profiles of a specified statistical quantity characterizing the GW field, such as horizontal velocity variance u 2 (e.g., Medvedev & Klaassen 1995), or vertical flux of horizontal momentum u w (e.g., Yigit et al. 2008). The former is convenient for comparison with observations of GW spectra. The latter is physically more lucid, because ρu w is an invariant in a nondissipative atmosphere. In GCMs, sources are specified by (1) prescribing the corresponding quantity at a certain level z s in the lower atmosphere, or (2) calculating it interactively using large-scale fields resolved by the model as an input. The latter is sometimes called "parameterization of gravity wave sources". Because mechanisms of wave excitation in the lower atmosphere are numerous, each requires a separate approach. To date, physically based schemes suitable for GCMs have been developed for GWs excited by convection (Chun & Baik 2002;Beres et al. 2004), flow over topography (McFarlane 1987), and fronts (Charron & Manzini 2002). In most other modeling studies, spectra at a source level are prescribed based on observational constraints, or simply tuned to obtain desired simulated fields. A comprehensive comparison of GW fluxes in observations and modeling has recently been performed by Geller et al. (2013). Although many GCMs use time-independent source spectra, GW excitation can undergo large changes during transient events, such as SSWs. Therefore, the importance of such variations should be explored and their possible impacts on the general circulation have to be taken into account in whole atmosphere GCMs.
In the middle atmosphere, the main mechanism of GW obliteration is nonlinear breaking and/or saturation that occurs when amplitudes become large. Therefore, most GW parameterizations developed for middle atmosphere GCMs (starting from that of Lindzen (1981)) have in common that they terminate harmonics, whose amplitudes reach a certain instability threshold. Exceptions are the approaches of Hines (1997) ("Doppler spread") and Medvedev & Klaassen (1995) ("nonlinear diffusion"), which sought to describe the underlying physics. The former is based on the assumption that harmonics are Doppler shifted by varying wave-induced wind directly to very short scales where they are removed by molecular diffusion. When averaged over wave phases, this parameterization, however, yields the very same termination of harmonics employing ad hoc chosen criteria. The approach of Medvedev & Klaassen (1995) is based on the concept of "enhanced diffusion" (Weinstock 1976;Weinstock et al. 2007). It takes into account Doppler shift by larger-scale harmonics in the spectrum, and erosion by shorter-scale ones. For parameterization purposes, Doppler shift can be neglected, the coefficient of eddy-induced diffusion is self-consistently cal-culated, and no "tuning parameters" are required (Medvedev & Klaassen 2000).
GW parameterizations suitable for thermosphere GCMs must account also for damping by molecular diffusion, thermal conduction, and ion friction. This is usually done by incorporating the respective dissipation terms into the complex dispersion relation in the form of imaginary parts of frequencies. The first parameterization of this kind has been proposed by Matsuno (1982), and the most recent derivation for molecular diffusion and thermal conduction has been performed by Vadas & Fritts (2005). This approach is based on the assumption that dissipation is relatively weak, where the degree of "weakness" depends on the characteristics of the harmonic and the background flow. This constitutes another limitation on GW parameterizations. Molecular viscosity grows exponentially with height in the thermosphere, and, eventually, the dissipation terms can significantly exceed all other balancing terms in the equations for waves. This means that GWs degenerate into other types ("viscous waves"), and can no longer be considered within the parameterization framework.
We illustrate the principles outlined above, and discuss some general details of implementation into a GCM using the extended GW parameterization (Yigit et al. 2008).
2.3. The extended nonlinear spectral gravity wave parameterization The word "extended" denotes that the parameterization has been extended to account for wave propagation in the thermosphere in accordance with the requirements outlined above (Yigit & Medvedev 2013). It solves the equation for the vertical structure of the horizontal momentum flux (per unit mass) u w associated with the harmonic j from a given spectrum of waves: where β j tot is the total vertical damping rate acting on the harmonic. If propagation is conservative (β j tot = 0), then the flux ρu w j is constant with height. The total damping rate for a given harmonic is the sum of the rates due to various dissipation processes affecting the propagation and acting simultaneously β j tot = β j non + β j mol + β j ion + β j rad + β j eddy + ... (2) The main processes accounted for by the scheme include, correspondingly, nonlinear breaking/saturation (β j non ), molecular diffusion and thermal conduction (β j mol ), ion friction (β j ion ), radiative damping (β j rad ), and eddy diffusion (β j eddy ) as suggested in the work by Yigit et al. (2008). The term β j non is parameterized after the work by Medvedev & Klaassen (2000), and comprises the effects of other harmonics on a given harmonic. Thus, the total wave field is not a simple collection of independent waves, but of interacting ones. The word "nonlinear" in the name of the parameterization signifies this. Dissipation of a harmonic is strongly affected by changes in the background wind as the vertical damping is inversely proportional to the intrinsic phase speed of the harmonic, i.e., β j ∝ (c j − u) −n , where the exponent n differs for various dissipation mechanisms (see e.g., Yigit et al. 2008Yigit et al. , 2012aYigit & Medvedev 2013). If the flux ρu w j changes with height, the wave momentum is trans-ferred to the mean flow by means of an acceleration or deceleration, which is often called "wave drag" The total "drag" is determined by the gradient of the sum of fluxes for all M harmonics in the spectrum, Σ M j a j . Equation (1) is solved for each grid column of a GCM. For that, values of u w j must be specified at a certain height z s in the lower atmosphere, which is considered as a source level. This initialization is done in all GW parameterizations, but the choice is extremely important for this scheme, because it contains no other tuning parameters, and the source spectrum is the only input. A representative spectrum can be seen in Yigit et al. (Figure 1, 2009), where the fluxes are specified as functions of horizontal phase velocities, and based on the observations of Hertzog et al. (2008). The "asymmetric" spectrum takes into account an anisotropy with respect to the mean wind at the source level. The latter has been first suggested heuristically (Medvedev et al. 1998), and a possible explanation has been offered recently (Kalisch et al. 2014).
GW harmonics with larger vertical wavelengths are less affected by dissipation and, therefore, tend to propagate higher. Typical scale height H also increases in the thermosphere (e.g., H is around 50 km at 250 km altitude). Because the parameterization is based on the WKB approximation (section 2.2), the vertical wavenumbers of accounted harmonics are limited by the relation k z H 1. This relation translates into the limitation on the maximum phase velocities of GW harmonics considered in the parameterization to be 80 to 100 m s −1 .
Using general circulation modeling, the extended GW scheme has been extensively validated against the empirical horizontal wind model (HWM) ) and the MSIS temperature distributions ). In a planetary atmospheres context, the extended scheme has successfully been used in a state-of-the art Martian GCM in order to investigate GW-induced dynamical and thermal coupling processes (Medvedev et al. 2013;Medvedev & Yigit 2012;Yigit et al. 2015;Medvedev et al. 2016).

EFFECTS OF INTERNAL GRAVITY WAVES ON THE
GENERAL CIRCULATION OF THE UPPER ATMOSPHERE Given the statistical approach to parameterizing waves, in which all the information on wave phases is lost, and given the set of assumptions listed in section 2.2, no effects of individual wave packets can be simulated with GCMs. They can only be approached with GW-resolving models similar to that of Miyoshi et al. (2014). Historically, the need for accounting for GW effects emerged from an inability of GCMs to reproduce the observed zonal mean circulation in the middle atmosphere (Holton 1983). In particular, the inclusion of parameterized effects of subgrid-scale waves has helped to realistically simulate the semi-annual oscillation in the MLT (mesosphere and lower thermosphere) with a GCM (Medvedev & Klaassen 2001). Manson et al. (2002) demonstrated the same for solar tides. Recently, Schirber et al. (2014) have shown that, with the use of a convection-based GW scheme, a GCM has reproduced a quasi-biennial oscillation (QBO) with realistic features.
Studying the effects of GWs of tropospheric origin in the thermosphere has a long history (see Yigit & Medvedev during high solar activity. For instance, the midlatitude easterly jet at 200 km is ∼20 m s −1 weaker in EXP2 than in EXP1, and the maximum of westerlies in the winter hemisphere of the same magnitude (∼30 m s −1 ) is shifted from 150 to 200 km.
[30] Distributions of the mean zonal torque by GWs are shown in Figures 10a and 10b. It is immediately seen that the main dynamical effects of GWs both in the middle and the upper atmosphere are to decelerate the mean zonal wind in an average sense. During the low solar activity, the GW drag is twice as strong in the high-latitude winter hemisphere and in the midlatitude summer one, and has approximately same magnitudes over the summer pole. The former result shows a good agreement with the findings of numerical studies conducted in the work by Vadas and Fritts [2006], who showed that the body force created from a localized convective plume has an amplitude twice as large during solar minimum as compared to solar maximum. The peaks in our simulations occur at different heights: around 200 and 280 km during small and large F 10.7 periods, respectively. GW momentum deposition is generally weaker below certain height during strong solar activity, which was also captured in our offline calculations. For instance, GW drag is weaker below ∼220 km at 70°S and 80°N during solar maximum, but exceeds the one for the solar minimum above this height. At 60°N, the wave-induced torque during solar maximum is greater only above 240 km, where all GWs from the source spectrum are almost entirely dissipated.
[31] The GW drag pattern looks somewhat differently when viewed in pressure coordinates (not shown here). Then, it gets shifted downward to higher pressures (lower logpressure heights) at high solar activity, and the magnitudes of the drag become weaker. This behavior is consistent with the notion that the molecular dissipation monotonically  wave numbers, respectively. Figure 3 shows the latitude-height distribution of the vertical convergence of u ′ w ′ . Below 100 km, the eastward (westward) wind acceleration occurs in the NH (SH), which acts to decelerate the stratomesospheric jets [Lindzen, 1981;Matsuno, 1982]. In the 100-120 km height region, the westward (eastward) acceleration appears in the NH (SH). The GW drag reaches À140 and 105 m s À1 (d) À1 in the NH and SH, respectively. The GW drag above 150 km height is much larger than that in the mesosphere and lower thermosphere (MLT). In particular, the eastward acceleration with a maximum value of +230 m s À1 (d) À1 is located at high latitudes. This result indicates that the GW drag plays an important role above the turbopause. Yiğit and Medvedev [2010] estimated the GW drag distribution in the thermosphere using CMAT2. Their result [Yiğit and Medvedev, 2010, Figure 10a] showed the strong eastward acceleration at middle and high latitudes between 150 and 250 km heights and is consistent with our result. However, our result is the first one obtained by a high-resolution GCM that can resolve GWs explicitly. The GW drag distribution in June solstice has been shown, and the distribution of the GW drag in December solstice is quite similar to that in June solstice (not shown).
The present GCM can resolve GWs with horizontal wavelength larger than 380 km. MF09 investigated the momentum flux u ′ w ′ distribution in the thermosphere as a function of the zonal wavelength. MF09 showed that the momentum flux due to GWs with wavelength of 500 km (850 km) is smaller than that of 1000 km by a factor of 3 (2). Thus, the momentum flux due to GWs decreases with decreasing wavelength of GWs. However, the momentum flux due to GWs with much shorter wavelengths is not well known and may affect the GW drag in the thermosphere. This means that the GW drag in the real atmosphere may be larger than the GW drag shown in Figure 3. In order to estimate effects of smaller-scale GWs on the GW drag, numerical simulation with higher horizontal resolution model is required.
The diurnal variations of the zonal wind and the GW drag are examined next. Figures 4a and 4b show the diurnal variations of the zonal wind at 20°S and 55°N, respectively. These representative latitudes are chosen because the diurnal mean GW drag in the thermosphere has maximum or minimum there. In this figure, the zonal wind is averaged for all longitudes, so that the zonal wind associated with the 2015, for more detail), however their dynamical importance at higher altitudes has not been fully recognized until recently. In all GCMs extending into the thermosphere, the effects of subgrid-scale GWs were either neglected, or assumed to decay exponentially above a certain height (e. g., turbopause ∼ 105 km). Simulations of  with the Coupled Middle Atmosphere and Thermosphere-2 ((CMAT2), Yigit 2009) GCM incorporating the extended nonlinear parameterization of Yigit et al. (2008) revealed that the momentum deposition by lower atmospheric GWs in the F region is substantial, and is comparable to that by ion drag. Figure 1 shows the latitude-altitude distribution of the simulated zonal mean zonal forcing by parameterized GWs. This forcing (known as "GW drag") is directed mainly against the mean zonal wind, and plays an important role in the momentum balance of the upper thermosphere, similar to the scenario in the middle atmosphere. The magnitude of thermospheric GW drag, exceeding ±200 m s −1 day −1 , is larger than its effects in the middle atmosphere. Miyoshi et al. (2014)'s recent simulations with a whole atmosphere GW-resolving GCM have confirmed 's predictions of the appreciable dynamical effects of lower atmospheric GWs on the general circulation of the thermosphere above the turbopause. Figure 2 presents the divergence of momentum fluxes (a in equation 3) due to the resolved portion of GW spectra (with horizontal scales longer than 380 km) calculated for solstice conditions (Miyoshi et al. 2014, Figure 3) as in the GCM modeling by . Considering the various approximations and limitations of the extended parameterization, and, especially, uncertainties with specifying GW sources, the two distributions in Figures 1 and  2 appear to be in a good qualitative and quantitative agreement. There are also some differences between the two simulations. In particular, in the Southern Hemisphere MLT, the high-resolution simulations show a region of eastward GW, which is only present at the Southern Hemisphere high-and low-latitudes in the parameterized simulation. Two possible sources of the discrepancies are the source spectrum and effects of the background winds on the propagation and the resulting dissipation. Overall, both simulation studies demonstrated that, due to propagation conditions in the middle atmosphere, most of the thermospheric GW activity concentrate at high-latitudes, where solar tides modulate local time variations of GW drag. This, and further analyses of the simulations with the high-resolution model provided evidences that thermospheric effects of GWs can be successfully parameterized in lower-resolution GCMs.
Thermal effects of GWs are two-fold: (a) heating due to conversion of the mechanical energy of dissipating harmonics into heat, and (b) heating and cooling associated with the downward sensible heat flux w T induced by these waves (Medvedev & Klaassen 2003;Becker 2004). Magnitudes of the former in the thermosphere are comparable with those due to the Joule heating, while the latter is comparable with the cooling rates due to molecular thermal conduction ), which suggests that the thermal effects of GWs cannot be neglected in the upper atmosphere. Yigit & Medvedev (2010)'s GCM simulations with the extended scheme have demonstrated that the variations of thermospheric GW effects are appreciable. GWs propagate to higher altitudes during high solar activity, but produce weaker drag than during periods of low solar activity. Their observations have later been qualitatively verified by the satellite observations of Park et al. (2014).

SUDDEN STRATOSPHERIC WARMINGS
4.1. Characteristics Sudden stratospheric warmings first discovered observationally by Scherhag (1952) are transient events during which the eastward zonal mean zonal winds weaken, or even reverse the direction at 60 • N (geographic) at ∼ 30 km (10 hPa), followed by the significant warming of the winter North Pole (90 • N) (Labitzke 1981;Andrews et al. 1987). Since the 1950s, as the interest in studying SSWs has grown, the classification of SSWs has evolved (see Butler et al. 2015, for a comprehensive discussion). Essentially, there are two commonly accepted types of warmings: a minor and a major warming. The warming is major if the equator-to-pole temperature gradient reverses poleward of 60 • latitude in addition to the reversal of the zonal mean zonal winds at 60 • N at 10 hPa (Labitzke 1981). If the westerly mean zonal wind weakens but does not reverse the direction, i.e., the stratospheric ionosphere [Forbes and Leveroni, 1992;Hagan et al., 2001;Abdu et al., 2006;Fuller-Rowell et al., 2008].
[4] A series of recent reports concerning ionospheric perturbations associated with lower atmospheric forcing have focused on sudden stratospheric warmings (SSW). These SSW events represent the most spectacular meteorological fluctuation in the polar stratosphere. During these events, stratospheric temperatures abruptly increase by tens of degrees, the normal winter polar vortex changes its position and shape, and zonal winds become weaker or even change direction [O'Neill, 2003]. Although not anticipated, observations have indicated that the ionospheric response to such events shows a global behavior. At middle latitudes, SSW is associated with warming in the lower thermosphere and cooling in the upper thermosphere, with both features exhibiting semidiurnal behavior [Goncharenko and Zhang, 2008]. At the magnetic equator, strong semidiurnal variation during the SSW event was observed in the vertical ion drifts [Chau et al., 2009]. This vertical motion leads to the large-scale redistribution of electron density in the daytime ionosphere, as demonstrated by Goncharenko et al. [2010] using total electron content (TEC) data measured by a global network of GPS receivers. An important aspect of this variation in TEC is its large magnitude, reaching 50%-150% compared to the pre-SSW behavior.
[5] In this paper, we analyze GPS TEC data for several recent stratospheric warming events, all of which occurred during time periods of low solar activity in the winters of 2007-2008 and 2008-2009. The key questions to be answered are the following: (1) whether the ionospheric variations reported by Chau et al. [2009] and Goncharenko et al. [2010] are repeatable for other stratospheric warming events; (2) what the characteristic signatures of ionospheric response to SSW in the low latitudes are; (3) what the temporal evolution of response to SSW in the low-latitude ionosphere is.

Data and Methods
[6] We use observations of global total electron content (TEC) obtained from the network of worldwide GPS receivers [Mannucci et al., 1998;Coster et al., 2003] that were calculated using the MIT Automated Processing of GPS (MAPGPS) software [Rideout and Coster, 2006]. The TEC estimates are produced in 1°× 1°bins of latitude/longitude with 5 min temporal resolution and distributed over those locations where data are available. The advantage to this process is that it is strictly data driven, with no underlying assumptions or models that smooth out real gradients in the TEC. The errors in the MAPGPS code are tracked throughout the processing, and random and correlated errors are handled separately. This allows optimal estimation of binned measurements using weighted averages and allows error values to be calculated independently for each binned measurement.
[7] This study includes GPS TEC data provided by the Low-latitude Ionospheric Sensor Network (LISN). LISN is a distributed observatory designed to monitor the state of the low-latitude ionosphere and to provide regional coverage in South America centered along the 70°W longitude. Recently installed LISN GPS receivers produce unprecedented spatial coverage in South America, allowing us to investigate the development of the equatorial ionization anomaly (EIA) in detail.
[    vortex does not break down, during a temperature increase at the Pole, then the warming is defined as a minor event.
An illustration of the major SSW features is seen in Figure 3 for a representative major warming that took place in the winter of 2008-2009, as adopted from the work by Goncharenko et al. (2010, Figure 1). These stratospheric conditions are based on data from the National Center for Environmental Prediction (NCEP). Within about 5 days, the zonal mean temperature at 10 hPa increases by more than 60 K (from 200 to more than 260 K) at the North Pole, that is, more than 30% increase (top panel). The average temperature at high-latitudes (60 • to 90 • N) increases significantly as well. The eastward (positive) zonal mean zonal wind starts decelerating already before the onset of the warming at the Pole and, reverses its direction, reaching a minimum over a period of about 10 days (bottom panel). The thin solid curves in each panel show the 30-year means of the associated parameters. Goncharenko et al. (2010) have also demonstrated in their analysis that the 2008-2009 warming was related to a weakening of the planetary wave-1 and an enhancement of the wave-2.
A comprehensive review of the earlier theoretical explanations of SSWs can be found in the works by Schoeberl (1978) and Holton (1980). Earlier studies have indicated that planetary-scale waves have to be properly taken into account during warming periods. According to the seminal work of Charney & Drazin (1961), planetary-scale disturbances can propagate from the troposphere into the stratosphere in the presence of prevailing westerlies, and the transport of eddy heat and momentum by vertically propagating waves is expected to modify the stratospheric zonal flow. Initial idealized simulations of wave propagation have suggested that planetary waves with wave numbers 1 and 2 can reach the stratosphere (Matsuno 1970). Matsuno (1971) modeled that Rossby wave-mean flow interactions decelerate the polar night jet, leading to weakening and even breakdown of the polar vortex, and ultimately to a sudden warming of the polar region. Later, the numerical works by Holton (1976) and Palmer (1981) have qualitatively provided supporting evidence for Matsuno (1971)'s model.

Mechanism of the sudden warming
In the winter (solstice) period, the Northern Hemisphere stratosphere is dominated by westerly jets whose strength increases with altitude. Quasi-stationary planetary waves can propagate vertically upward, provided that the mean zonal flow satisfies the conditions for vertically propagating wave modes. For these waves, the zonal wind has to fulfil the following condition (Holton & Hakim 2012, Equation (12.16)): where the Rossby critical velocity u c is defined in terms of the characteristics of the background atmosphere and wave by where k 2 h = k 2 + l 2 is the horizontal wavenumber that depends on the zonal (k = 2π/λ x ) and the meridional (l = 2π/λ y ) wavenumbers; f = f 0 − βy is the beta-plane approximation for the Coriolis parameter, and β ≡ ∂f ∂y is the meridional gradient of the Coriolis parameter. The condition (4) suggests that planetary waves can propagate vertically only in the presence of westerly winds that are weaker than a certain critical value u c , which depends on the horizontal scale of the wave. Dynamical conditions are, therefore, favorable for the vertical propagation of planetary waves in the winter Northern Hemisphere with prevalent mean westerly winds. This condition is important for understanding the propagation of GWs, which are also affected by the mean wind distributions. Namely, before the warming, the stratospheric zonal mean winds are eastward. They filter out a significant portion of the eastward directed GWs, favoring the upward propagation of harmonics with phase velocities directed westward. During the warmings, the decelerating westerlies increase the chances of GWs with eastward horizontal phase speeds to propagate to higher altitudes (Yigit & Medvedev 2012).
In the winter stratosphere, waves are rapidly attenuated, thus decelerating the mean zonal flow. For the occurrence of SSWs, a large-scale wave transience, in particular, rapid temporal changes of planetary wave activity are also important. They maintain the convergence of the westward momentum flux, leading to strong polar night jet deceleration and poleward meridional flow enhancement (Andrews et al. 1987). Additionally, radiative forcing sustains a cold winter North Pole with negative equator-to-pole mean temperature gradient, that is, ∂T ∂y < 0. The rapid deceleration of the stratospheric mean flow implies a decreasing (positive) vertical gradient of the zonal flow between the troposphere and stratosphere. From the thermal wind relation ∂u ∂z ∼ − ∂T ∂y , this decrease implies a rise of temperature at the winter pole, meaning that the equator-to-pole mean temperature gradient becomes less negative. During a major warming, this gradient even reverses due to the reversal of the vertical gradient of zonal mean wind. The strong polar night jet deceleration leads to a departure from the thermal wind balance, and the poleward meridional flow, which is caused by the Coriolis force associated with the westward forcing, is induced to recover this balance. This enhancement of the Brewer-Dobson circulation ultimately results in an adiabatic warming at Northern Hemisphere high-latitudes.

OBSERVED CHANGES IN THE UPPER ATMOSPHERE DURING SUDDEN STRATOSPHERIC
WARMINGS Given the rapid and strong local changes in the circulation and thermal structure of the stratosphere during SSWs, the natural questions that bear in mind are (1) how high the effects of the warming propagate in altitude, and (2) to what extent the changes in the upper atmosphere can be associated with the sudden warmings. Planetary waves cannot propagate directly to much higher altitudes, but the stratosphere and mesosphere are closely connected via circulation and by GWs and tides. As sudden warmings and the associated dynamical changes in the stratosphere occur over relatively long time scales (e.g., ∼ 10 days) compared to the periods of internal waves, lower atmospheric wave disturbances have sufficient time to propagate to higher altitudes, provided that propagation conditions are favorable. Therefore, one ought to expect a certain degree of coupling between the stratosphere and higher altitudes, probably beyond the middle atmosphere.
How can one associate observed upper atmospheric changes with SSWs? Essentially, a ground-to-upper atmosphere observation with a single instrument is beyond the capabilities of the current technology. For the purposes of observational analysis, SSW events/periods ought to be identified. For this, an appropriate description of stratospheric dynamics is needed in the first place. This representation could be, for example, obtained from numerical forecast models that assimilate in-situ and remote-sensing data, such as the European Centre for Medium-Range Weather Forecast (ECMWF) analyses, and produce the required global fields of atmospheric parameters. Then, observational data can be investigated together with the numerical model output (e.g., Pancheva et al. 2008).
The deceleration of the stratospheric eastward zonal flow during sudden warmings leads, ultimately, to an upward circulation in the mesosphere that results in mesospheric cooling (Liu & Roble 2002). Such direct link between these two regions have motivated a number of scientists to investigate the details of stratosphere-mesosphere changes during warmings. Based on temperature and geopotential height data obtained from the Sounding the Atmosphere using Broadband Emission Radiometry (SABER) instrument of the Thermosphere Ionosphere Mesosphere Energetics and Dynamics (TIMED) satellite and the VHF radar horizontal winds, Pancheva et al. (2008) have investigated planetary wave-induced coupling in the stratosphere-mesosphere during the major warming of 2003/2004 winter Northern Hemisphere. Yuan et al. (2012) studied the response of the middlelatitude mesopause region to the 2009 major SSW, using a sodium Doppler wind-temperature lidar. They have discovered anomalous behavior of the mean temperature and zonal winds around the mesopause during the warming, and have concluded that it was due to a direct impact of the major warming on the middle-latitude mesopause. The 2009 SSW has been one of the strongest warming events that has been recorded. The features around the mesopause during SSWs can be largely characterized in terms of an "elevated stratopause", which forms around 75-80 km after the SSW occurrence and then descends (Manney et al. 2009). The role of GWs and planetary-scale waves in the time evolution of the elevated stratopause have been investigated by a number of authors (e.g., Siskind et al. 2010;Chandran et al. 2011;Limpasuvan et al. 2012).
Vertical coupling between the stratosphere and the lower thermosphere has been studied in the low-and middle-latitude Northern Hemisphere winter of 2003/2004 based on the temperature data from SABER/TIMED (Pancheva et al. 2009). According to Goncharenko & Zhang (2008)'s analysis of the Millstone Hill incoherent scatter radar (ISR) ion temperatures data, warming in the lower thermosphere and cooling above 150 km were revealed during a minor SSW. Using data from the Jicamarca ISR, Chau et al. (2009)  Based on TEC (total electron content) data retrieved from a worldwide network of GPS observations, Goncharenko et al. (2010) have found a significant local time modulation of the equatorial ionization anomaly (EIA) induced by SSWs. Using the European Incoherent Scatter (EISCAT) UHF radar, Kurihara et al. (2010) have detected short-term variations in the upper atmosphere during the 2009 major warming. In their analysis of Fabry-Perot and incoherent scatter radar data, Conde & Nicolls (2010) have identified that the period of reduced neutral temperatures at 240 km, which corresponded closely to the main phase of the warming.
More recently, analyzing the Constellation Observing System for Meteorology Ionosphere and Climate (COSMIC) data, Pancheva & Mukhtarov (2011) have found a systematic negative response of ionospheric plasma parameters (f 0 F 2 , h m F 2 , and n e ) during the warmings of 2007-2008 and 2008-2009. An illustration of their results for the mean zonal mean electron density (in MHz) at 300 km are seen in Figure 4, where the 2007-2008 and 2008-2009 warming events are shown on the left and right panels, respectively. The response to the warming is negative and mainly occurs in the low-and middle-latitude region. Liu et al. (2011) used neutral mass density observations from the CHAMP (Challenging Minisatellite Payload) and GRACE (Gravity Recovery and Climate Experiment) satellites to study the thermospheric variations during the 2009 major warming. They have found a substantial decrease of the mass density, and concluded that it was potentially associated with a strong thermospheric cooling of about 50 K. Goncharenko et al. (2013) have investigated the day-to-day variability in the middle-latitude ionosphere during the 2010 major warming using the Millstone Hill ISR. They have discovered that semidiurnal and terdiurnal tidal variations were enhanced during the SSW. Jonah et al. (2014) have used a suit of observational data from GPS, magnetosphere, and meteor radar in order to investigate the response of the magnetosphere and ionosphere to the 2013 major SSW. Analyzing long-term data of the global average thermospheric total mass density derived from satellite orbital drag, Yamazaki et al. (2015) showed density reduction of 37 % at 250575 km during SSW period that can be associated with lower atmospheric forcing. Recently, using data from GPS and ionosonde stations, Fagundes et al. (2015) investigated the response of the low-and middle-latitude ionosphere in the Southern Hemisphere to the 2009 major SSW and found that during the warming, TEC was depressed following the SSW temperature peak.
Overall, these studies suggest that (1) a variety of instru-(left column of plots) and 2008/2009 (right column of plots). Again the ionospheric response is different for both winters. While the response in 2007/2008 is composed of several decreases of the electron density centered near days 60, 90, 115 and 140 (responses to the temperature peaks from Fig. 1a-left plot), the winter 2008/2009 consisted of mainly one electron density reduction centered near day 120, which is a few days delayed response to the temperature peak from Fig. 1a (right plot). Fig. 2a shows two main features of the latitude structure of the mean electron density response: (i) the ionospheric response amplifies toward the equator and (ii) it is not symmetrical with respect to the modip latitude; the response is mainly in the Northern Hemisphere (NH).  Fig. 1a. The altitude structure of the mean electron density shows that the response is observed mainly above the F-region maximum (above $ 300 km height).
We mentioned before that the DW1 electron density oscillation is produced by a combination of the diurnal variability of the photoionization and the effect of the migrating diurnal tide forced from below and that the used data analysis method cannot separate the two effects. However, Mukhtarov and Pancheva (2011) presented evidence indicating that the DW1 electron density oscillation is shaped mainly by the photoionization, i.e. is due to a solar zenith angle effect. This is the reason for this oscillation to be called diurnal variability of the ionosphere. Fig. 3a shows the latitude structure of the COSMIC electron density diurnal variability at an altitude of 400 km in winters 2007/2008 (left plot) and 2008/2009 (right plot). Again different temporal variabilities of the ionospheric response in both winters is detected. Several falls of the diurnal amplitude in the first winter, centered near days 60, 90, 120 and 145 (very weak), can be seen. A single strong reduction of the diurnal amplitude in the second winter, centered near day 120, can be distinguished; it is preceded by an increase of the diurnal amplitude centered near day 105. The average reduction of the electron density at 400 km height during both winters is $0.3 MHz, which means $11% change of the electron density. The latitude structure of the COSMIC diurnal variability clearly indicates a response only in the NH. . The temporal variability of the diurnal response shows the same feature as its latitude structure. The altitude structure of the diurnal response, particularly that in the winter of 2008/2009, indicates that the response happens mainly above the F-region electron density maximum, e.g. above 300-350 km height.

Discussion and conclusions
In this paper we have presented the global spatial (latitude and altitude) structure and temporal variability of the mean ionospheric response to SSW events in winters of 2007/2008 and 2008/2009. To elucidate the effect of the SSWs on the ionosphere the COSMIC f o F2, h m F2 and electron density data at fixed altitudes are analyzed in terms of their mean (zonal and time) values, various migrating and nonmigrating tidal and SPW components (see formula (1)). The results only for the zonal and time means (Y 0 from formula (1)) and diurnal (migrating) variability have been reported here. As presented in Fig. 1b, the negative response to the SSW events of the mean f o F2 and h m F2 parameters indicates that such ionospheric variability can be detected by ionosonde measurements as well. Many years ago Pancheva and Spassov (1989) investigated the effect of the SSWs on the ionospheric D-and F-region variabilities by analyzing the radio-wave absorption and ionosonde f o F2 data, respectively. The authors found negative responses of both regions; usually the F-region response lags $ 1-2 days behind that of the D-region.    ments has been used to study upper atmospheric changes during SSWs; (2) They convincingly demonstrated that SSW events affect the thermosphere-ionosphere system beyond the turbopause; and (3) the associated observed changes in the upper atmosphere vary from one warming event to another. Some studies indicate that large-scale internal wave processes may be involved in vertical coupling during SSWs. One of the less appreciated topics in SSW studies is the role of smallscale GWs. We next discuss the upper atmosphere changes during SSWs in the context of lower atmospheric small-scale GWs that can propagate to the thermosphere (Yigit et al. , 2012a. 6. UPPER ATMOSPHERIC CHANGES DURING SUDDEN STRATOSPHERIC WARMINGS Observing dynamical changes, e.g., with satellites and radars, cannot provide sufficiently detailed information on characteristics and physical mechanisms of vertical coupling. Observations may, and in fact do raise new questions, which can be investigated by models. A powerful tool is general circulation models (GCMs) that solve the coupled governing equations of atmospheric and ionospheric physics in time and three-dimensional space. GCMs generate a full set of field parameters that can be diagnosed in detail in order to investigate the physical mechanisms that shape the state and evolution of the atmosphere. Therefore, global models can provide an unprecedented insight in vertical coupling processes between the different atmospheric regions. One should nevertheless be aware of the limitations of GCMs, such as, resolution, boundary conditions, and parameterizations.
As discussed in section 3, GWs have a profound effect on the dynamical (Yigit et al. , 2012aMiyoshi et al. 2014;Vadas et al. 2014), thermal Hickey et al. 2011) and compositional (Walterscheid & Hickey 2012) structure of the upper atmosphere. The state of the back-ground middle atmosphere plays a crucial role in modulating the propagation of GWs into the thermosphere. Given that SSWs modify the middle atmospheric circulation significantly, how can they influence the upper atmosphere in the context of GW propagation and dissipation in the thermosphere? Resolving this science question requires a use of comprehensive GCMs with appropriate representation of small-scale GWs. The general circulation modeling of Yigit & Medvedev (2012) using the extended GW parameterization of Yigit et al. (2008) has demonstrated GW propagation and appreciable dynamical effects in the upper thermosphere during a minor warming. The universal time (UT) variations of the GW-induced zonal mean root-mean-square (RMS) wind fluctuations (in m s −1 ) and zonal mean GW drag (in m s −1 day −1 ) are shown in Figures 5a,b. GW-induced RMS wind fluctuations are given by where M is the number of harmonics in the spectrum and variance u 2 j is related to the wave amplitude as u 2 j ≡ |u j |. The GW RMS wind fluctuations are a proxi for the subgridscale GW activity as the fluctuations induced by all waves in a GW spectrum are taken into account and do not represent the resolved wind fluctuations. In the course of the warming, GW activity increases by a factor of three exceeding 6 m s −1 in response to weakening of the polar night jet. In addition to persistently large values in the lower thermosphere, modulation of the GW activity is seen higher in the thermosphere. Following the increase of GW activity, (eastward) GW drag increases in the thermosphere during the warming as well.
The effects of GWs in the upper atmosphere during SSWs are not confined to only those in a zonal mean sense. Recently, Yigit et al. (2014) have investigated the details of GW temporal variations in the thermosphere during a minor warming simulated with a GCM. They modeled that GW drag and SSW, as defined in the seminal review of Schoeberl [1978]: "In a minor warming, …the meridional temperature gradients usually weaken but do not reverse." In our simulations, T increases rapidly by 18 K from 216 K on Dec 16 to 234 K on 4 Jan, and steadily decreases again after this date. D T decreases throughout the period when T at the North Pole rapidly rises. This reduction of the equatorward temperature gradient occurs due to the poleward transport of heat by enhanced planetary wave activity. Figure 1b shows that before the minor warming commences, zonal mean zonal wind u weakens, and decreases continuously from 52 to 24 m s À1 throughout the minor warming.
[14] To demonstrate the spatial evolution of the fields during the warming, the NH polar stereographic projections of geopotential height (color shaded) and zonal wind (contour) at 10 hPa are shown on three dates (15 Dec, 24 Dec, and 1 Jan) in row two of Figure 1. As the warming progresses, the polar night vortex splits into two and the geopotential height increases at high-latitudes in accordance with the rising temperature.

Variations of Gravity Wave Activity and Effects in the Thermosphere
[15] Next, we investigate how the minor SSW shown in Figure 1 impacts the GW propagation/dissipation, and what effects these waves produce above the turbopause ($105 km). We evaluate the root mean square (RMS) wind , M being the number of harmonics in the spectrum, s 2 is the horizontal wind variance, and u′ i is the wave amplitude in the zonal wind, as a quantitative measure of wave penetration into the upper atmosphere. GW dynamical effects are characterized by the zonal GW drag, a x ¼ r À1 d ru′w′ À Á =dz, where r is the mean density.
[16] Figure 2a presents temporal variations of the zonally averaged RMS at 60 N from 14 Dec to 7 Jan as a function of altitude. The white dashed line marks the onset of the warming when the westerly zonal wind started to weaken (cf. Figure 1b). GW-induced wind variations above 90 km steadily rise after that date, and reach their maxima around 25 Dec. The white dotted line on Jan 4 indicates the peak of the SSW. Before the warming, the maximum of RMS ($3 m s À1 ) is in the mesosphere at around 80 km. During the warming, it increases to $6 m s À1 , and is shifted higher into the lower thermosphere to around 120 km. Also, the overall GW penetration into the thermosphere increases during the event. Thus, the RMS rises to > 4 m s À1 at altitudes up to $250 km compared to 1-1.5 m s À1 before the onset.
[17] UT variations of the mean zonal GW drag a x at 60 N are plotted in Figure 2b. Remarkable changes occur with it during the SSW. Penetration higher into the thermosphere and larger amplitudes of waves result in the enhancement of a x , especially in the lower thermosphere at around 120 km, and in the upper thermosphere at $250 km. Before the SSW, there is only easterly drag below 120 km with the maximum of 100 m s À1 day À1 at around 80 km, and a weak (up to 10 m s À1 day À1 ) westerly drag in the upper thermosphere at 60 N. Mean zonal wind u shown in Figure 2c demonstrates the role of selective filtering in GW effects in the thermosphere. After the warming onset, stratospheric westerlies weaken, thus permitting more eastward harmonics to propagate upward. As a result, the deposited momentum is markedly westerly between 120 and 180 km. Higher in the thermosphere at around 250 km, the westerly drag dramatically enhances after 21 Dec to more than 150 m s À1 day À1 .  its variability is dramatically enhanced in the thermosphere during the warming and thus lead to a ∼ ±50% modulation of small-scale and short-term variability in the resolved thermospheric winds, where the small-scale variability has been evaluated by subtracting the contributions of the large-scale tides. Recently, Miyoshi et al. (2015) have demonstrated with a GW-resolving GCM that the SSW has major dynamical and thermal impact on the upper atmosphere, substantially influencing the global circulation. Changes in the mean zonal wind produce a feedback on GWs by modifying filtering, dissipation, and propagation conditions. The upper atmosphere above the turbopause has a great amount of variability owing to the simultaneous influences of meteorological and space weather processes (Matsuo et al. 2003;Anderson et al. 2011;Yigit & Ridley 2011b;Yigit et al. 2012b). Often, separating the components and sources of variability in observations is a challenging task. Thus, following their observations of an SSW, Kurihara et al. (2010) have concluded that understanding the link between SSWs and thermal and dynamical changes in the upper atmosphereionosphere requires investigations of GW-mean flow interactions processes. GCM studies can greatly supplement these efforts.
Predictions with GCMs indicate that small-scale GWs can substantially contribute to the variability of the upper atmosphere. Also, recent modeling studies with a wholeatmosphere GCM have shown an enhancement of the semidiurnal variation in the ionospheric E×B drifts during the 2009 major warming (Jin et al. 2012). This increase has been interpreted as a consequence of the semidiurnal tidal amplification in the lower atmosphere. Further investigations that incorporate a fully two-way coupled thermosphere-ionosphere under the influence of lower atmospheric waves are required in order to assess the significance of the lower atmosphere in the context of upper atmosphere variability. In characterizing the upper atmosphere processes, the variability is always defined with respect to some appropriate mean. Therefore, the quantity of variability is not uniquely defined, and care should be taken when comparing one study to another. In a broader context, the presence of any kind of variability restricts scientists' ability to predict the future state of the atmosphere, and it is crucial to determine the sources of variability and quantify the magnitude thereof. 7. CONCLUSIONS This paper has briefly reviewed the current state of knowledge and most recent developments with understanding the role of GWs in vertical coupling during SSWs. The observed upper atmosphere changes during SSWs have been presented. An emphasis was placed on the processes above the mesopause, and on how they can be studied with general circulation models.
The geosciences community increasingly recognizes that the effects of lower atmospheric gravity waves extend beyond the middle atmosphere into the atmosphere-ionosphere system and are of global nature. Similarly, sudden stratospheric warmings were used to be looked upon as stratospheric phenomena, but now compelling observational evidences of their signatures in the thermosphere-ionosphere are being routinely provided.
With the rapid progress in the field of atmospheric coupling, further key science questions on the role of GWs in coupling atmospheric layers arise: • What are the spectra of gravity waves in the lower and upper atmosphere? How do they change during the different phases of SSWs?
• How well do GW parameterizations describe wave spectra and reproduce their effects during SSWs?
• What is the relative role of GW-and electrodynamical coupling between atmospheric layers in the variability of the atmosphere-ionosphere system during SSWs?
• What are the effects of GWs on the circulation and thermal budget of the upper atmosphere during major sudden stratospheric warmings?
• Do GWs in the upper atmosphere affect the development of sudden stratospheric warmings, or they are a mere reflection of processes in the lower atmosphere?
• Do GWs have a role in latitudinal coupling in the thermosphere during SSW events?
This is certainly an incomplete list of scientific questions, answering which requires concerted observational, theoretical, and modeling efforts from scientists of both lower and upper atmosphere communities.