Gravity Waves induced Wind Shears Derived from SABER Temperature Observations

. Large wind shears around the mesopause region play important roles in atmospheric neutral dynamics and ionospheric electrodynamics. Based on previous observations using sounding rockets, lidars, radars and model simulations, 15 large shears are mainly attributed to gravity waves (GWs) and modulated by tides (Liu, 2017). Based on the dispersion and polarization relations of linear GWs and the SABER temperature data from 2002 to 2019, a method of deriving GW-induced wind shears is proposed. The zonal mean GW-induced shears have peaks (13-17 ms -1 km -1 ) at around the mesopause region, i.e., at z=90-100 km at most latitudes and at z=80-90 km around the cold summer mesopause. This latitude-height pattern is robust over the 18 years and agrees with model simulations. The magnitudes of the GW-induced shears exhibit year-to-year 20 variations and agree with the lidar and sounding rocket observations on climatology sense but are 60-70% of the model results in the zonal mean sense. The GW-induced shears are hemispheric asymmetric and have strong annual oscillation (AO) at around 80 km (above 92 km) at the northern (southern) middle and high latitudes. At middle to high latitudes, the peaks mean shears. The maxima of the std and the top 10% largest shears are, respectively, ~12 ms -https://doi.org/10.5194/acp-2020-515Preprint.Discussionstarted:1July2020c The GW-induced shears derived here agree with previous lidar and sounding rocket observations in the aspects of height structures and magnitudes on the climatology sense. Moreover, the GW-induced shears derived here agree with the high-resolution model simulation results in the aspects of latitude-height patterns but have smaller magnitudes in the zonal mean sense. The GW-induced shears reach their maxima around the mesopause region and increase with the increasing latitudes. At most latitudes and during all seasons, the maxima of GW-induced shears are at ~ z =90-100 km. At high latitudes of the summer 5 hemisphere, the maxima of GW-induced shears at a lower height (~ z =80-90 km). This latitude-height pattern of GW-induced shears is independent on years. The magnitudes of the GW-induced shears exhibit year-to-year variations. The GW-induced shears exhibit more prominent AO and SAO at high latitudes than those at lower latitudes. The height variations of the amplitudes AO and SAO are hemispheric asymmetric. The strong AO occurs at around 80 km in the NH and above 92 km in the SH. At middle to high latitudes, the phases of AO shift from winter to summer and then to winter again 10 with the increasing height. The amplitudes of SAO decrease with the decreasing latitudes. The phases of SAO are in May and June when the SAO reach its peak at middle to high latitudes. The main limitation of the method is the overestimation of the GW-induced shears due to the unresolved GW propagation direction by the method. The other limitations, such as the observational filter, long sampling distance and cutoff criterion of the vertical wavelength, will underestimate the GW-induced shears. To overcome these limitations, it is necessary to develop new techniques and/or more dense satellites to observe the atmosphere with higher horizontal and vertical resolutions.

the ALOMAR observatory at 69.3°N) in July 2002, Fritts et al. (2004) showed large winds and shears at z~85-95 km and ascribed them to strong gravity waves (GWs) activity. These GWs were generated by convection, orography, jet and front in the troposphere and stratosphere and had increasing wave amplitude and shears when they propagated into the lower thermosphere. It is now accepted that the large winds and shears are a common phenomenon in the MLT region (Liu, 2007;2017). The large winds and shears play important roles in forming the middle latitudes sporadic E layers and driving the 5 equatorial electrojet (Mathews 1998;Hysell et al., 2002;Haldoupis, 2012;Shinagawa et al., 2017;Yu et al., 2019), in controlling atmospheric stabilities and the propagation of GWs, and in transporting and mixing tracers locally and/or globally (Fritts et al., 2004;Liu 2007Liu , 2017Yue et al., 2013;Stevens et al., 2014).
Based on the theory of dynamic instability, Liu (2007) showed that the maximum wind shears derived from = 2 ( : the buoyancy frequency, : the vertical shear) agree well with the observed large wind shears. However, the global scale 10 models (e.g., TIME-GCM: Thermosphere, Ionosphere, Mesosphere, Electrodynamics General Circulation Model) cannot reproduce the observed large winds and wind shears but can increase the amplitudes of winds and wind shears with finer spatial resolutions (Larsen & Fesen, 2009). This indicates that the model resolution and thus tides and small-scale waves (e.g., GWs) are important in driving large winds and shears. Using a local numerical model, Liu et al. (2014a) showed that the nonlinear interactions between GWs and tides can produce large winds and shears in the MLT region. The tidal phases modulate the 15 peak height of large winds and shears. Using the spectral element version of NCAR Whole Atmosphere Community Climate Model (WACCM) with horizontal resolution of ~25 km and vertical resolution of 0.1 scale height, Liu (2017) reproduced the large wind shears, which are in good agreement with observations. Through scale separation, Liu (2017) proposed that smallscale waves (with zonal wavenumber > 6), likely GWs, play a dominant role in producing large shears. The high resolution WACCM can resolve GWs with scales longer than ~100 km. Tidal waves make secondary contribution to the magnitudes of 20 shears but can modulate the shears produced by GWs.
Large winds and shears have been observed both by sounding rocket and ground based lidar and radar, all locally and cannot provide a global morphology. High resolution GCM (e.g., WACCM) simulations can provide a global picture of large winds and shears but need be validated through global observations (e.g., satellites). Moreover, it is still a challenging to study the intra-annual and/or the inter-annual variations of large winds and shears through high resolution GCM simulations due to 25 the computational cost (Liu, 2017(Liu, , 2019.
Satellite observations provide a good opportunity to study the climatology of global winds. The Wind Imaging

Interferometer (WINDII) and the High Resolution Doppler Imager (HRDI) instruments onboard the Upper Atmospheric
Research Satellite (UARS) provide global observations of winds in the MLT region (McLandress et al, 1996;Zhang et al., 2007;Shepherd et al., 2012). Combining the winds observed by HRDI and data assimilation system, Swinbank & Ortland 30 (2003) developed a climatology that describes the monthly zonal mean zonal winds from the surface to the upper mesosphere.
The TIMED Doppler Interferometer (TIDI) instrument on board the Thermosphere Ionosphere Mesosphere Energetics and Dynamics (TIMED) satellite measures global winds in the MLT region Niciejewski et al., 2006). These observations advanced our knowledge on the global winds in the MLT region. However, the global characteristics of shears https://doi.org /10.5194/acp-2020-515 Preprint. Discussion started: 1 July 2020 c Author(s) 2020. CC BY 4.0 License. are poorly known. Since large shears play important roles in the atmospheric dynamics and in the MLT region and ionospheric E region and are likely caused by GWs (Fritts, 2004;Liu, 2007Liu, , 2017Yue et al., 2010), it should be possible to derive wind shears by combining the GW theory and GWs derived from other observed physical quantities (e.g., temperature).
The remainder of this paper is organized as follows. Section 2 presents the method of deriving wind shears (short for 10 GW-induced shears) induced by GWs, and the uncertainties will be presented. Section 3 presents the comparisons of the GWinduced shears with the model and observational results. Then the latitudinal and intra-annual variations of the GW-induced shears are presented in Section 4. The limitations of the method and their possible influences on the GW-induced shears are discussed in Section 5. A short summary is given in Section 6.

Theory of Deriving GW-induced Shears
The basic idea of deriving GW-induced shears is the linear GW theory, which includes the dispersion and polarization relations of a monochromatic GW. A monochromatic linear GW has the form of (Fritts & Alexander, 2003), here, = √ −1 is the imaginary unit. The subscript denotes a monochromatic GW, ′ and ′ are the horizontal wind 20 perturbations. ′ and � are the perturbation and background temfperatures, respectively. � , � and � are the amplitudes of ′ , ′ and ′ � ⁄ , respectively. = + + − is the phase of GW. , and are the wavenumbers in the horizontal ( , ) and vertical (z) directions, respectively. and t are the ground-based frequency and time.
Based on the polarization of the monochromatic linear GWs with lower and medium frequencies (Fritts & Rastogi, 1985;Fritts & Alexander, 2003), the relations between � , � and � can be derived as (Eckermann et al., 1995;Liou et al., 2003;25 Gubenko et al., 2008), here � and = 2Ω sin (Ω = 7.292 × 10 −5 s −1 , is latitude) are the intrinsic and inertial frequencies, respectively. The wind shear of each monochromatic GW can be written as, In real atmosphere, a GW profile consists multiple spectra and can be regarded as a superposition of monochromatic GWs.
For each monochromatic GW, Eq. (2-5) are valid. Here, we use ′ to represent observed GWs, which contains multiple monochromatic GWs (e.g., many ′ ) and can be expressed as, 5 In the same way, the component ′ of an observed GWs can be expressed as, For the GW-induced shears of ′ and ′ , we have Finally, the GW-induced shear can be written as

GW-induced Shear Derived from Synthesized GW Profiles and Validations
To demonstrate the applicability of the theory and the procedure of retrieval the GW-induced shears, we construct two 15 synthesized temperature perturbation profiles (shown in Fig. 1a). Each profile is the sum of three monochromatic GWs with three different vertical wavelengths and height dependent amplitudes, The background temperature � ( ) = 240 K. The black line in Fig. 1a show the 1 ′ ( ) and noted as 1 ′ . The profile 2 ′ (red line in Fig. 1a) has the same amplitudes and vertical wavenumbers as 1 ′ but setting 1 = 2 = 3 = 2 ⁄ . 1 ′ ( ) and 20 2 ′ ( ) are used to represent two adjacent SABER measurements. It should be noted that we set the three monochromatic GWs in 2 ′ ( ) having the same phase shift of 2 ⁄ only for the convenience of theoretical representation. In real atmosphere, GWs with different vertical wavelengths may have different horizontal wavelengths and thus different phase shifts. These phase shifts can be calculated by comparing the phases of the two monochromatic GWs with same vertical wavelength embedded in the two adjacent GW profiles. While the phase of each monochromatic GW can be derived through discrete Fourier 25 transformation (DFT, Press et al., 1992). Then the phase shifts can be used to estimate horizontal wavenumber (e.g., Preusse https://doi.org/10.5194/acp-2020-515 Preprint. Discussion started: 1 July 2020 c Author(s) 2020. CC BY 4.0 License. et al., 2002;Ern et al., 2004;Alexander et al., 2008Alexander et al., , 2018Wang and Alexander, 2010). For example, if we assume the horizontal distance (∆ ) between the two profiles is 300 km, the phase shift of ∆ = 2 ⁄ indicates the horizontal wavenumber For the lower and medium frequency GWs, the dispersion relation can be simplified as (Fritts & Alexander, 2003), After we get the horizontal wavenumber ,ℎ from a GW profiles pair, then the intrinsic frequency for ′ ( ) can be calculated by Eq. (12). Using Eq. (2-5) and the prescribed amplitude and vertical wavenumber of each monochromatic GW, we can get the zonal and meridional winds components and their shears. Then, according to Eq. (6-9), the vector sum of the three monochromatic GW-induced wind profiles can be obtained and are shown as black lines in Fig. 1b for zonal wind and 1(c) for meridional wind. Here we assume the latitude is at 30ºN, which is a typical latitude at middle latitudes. The corresponding 10 shears are also shown as black lines in Fig. 1d-e. Since the amplitude, vertical and horizontal wavenumbers of GW profile pairs are prescribed and are not derived through DFT, which is a key step of the spectral decomposition method described below, the retrieved winds and shears are referred to as theoretical results. These theoretical results can be used to measure the results obtained from the spectral decomposition method proposed below. It should be noted that the zonal and meridional winds are the winds along and cross the orbit track directions, respectively, and may not coincide with the eastward and 15 northward directions, but are the winds along and cross the orbit track directions, respectively. Here zonal direction is points from the location of 1 ′ ( ) to that of 2 ′ ( ). The meridional direction is 90° counterclockwise from the zonal direction.
Now we describe the method of retrieving winds and wind shears induced by GW profile pairs, whose amplitudes, vertical and horizontal wavenumbers are not prescribed but should be evaluated. We named this method as "spectral decomposition method" since the principle ideas are: (1) decomposing an observed GW profile into multiple monochromatic waves; (2) 20 applying linear GW theory on each monochromatic wave to get monochromatic wind and shear of each wave component; (3) finding the vector sum of monochromatic winds and shears to get the wind perturbations and shears induced by the observed GW. The detailed application of this method is described as the following three steps.
The first step is to evaluate the amplitude and vertical wavenumbers of each GW profile by the method of DFT (Press et al., 1992). The amplitude ( � , 1 ) and phase ( , 1 ) at the vertical wavenumber ( ) can be calculated through, 25 here the vertical wavenumbers = 2 (90 km) ⁄ , which corresponds the vertical wavelength ranging from ~5 km to ~30 km for a vertical extent of 90 km, which is the height coverage (18-108 km) of the SABER temperature profiles. In a same manner, we can get the amplitude ( � , 2 ) and phase ( , 2 ) of GW profile 2 ′ ( ) at .
The second step is to evaluate the horizontal wavenumber through the phase shift between two adjacent GW profiles. For 30 each Fourier component of , the phase shift is ∆ = , 1 − , 2 . According to the distance between the two adjacent profiles ∆ = 300 km, we can get ,ℎ = ∆ ∆ ⁄ = �∆ � (300 km) ⁄ . Then the intrinsic frequency � for the component of can be calculated by Eq. (12). Here we note that the horizontal wavelengths of GWs in real atmosphere may be shorter than 2∆ , only GWs with horizontal wavelengths longer than 2∆ are considered due to the sampling distances and the limb scanning mode of the SABER instrument (Preusse et al., 2002;Ern et al., 2004).
The third step is to calculate the GW-induced wind (shown as red dashed lines in Fig. 1b-c) and shears perturbations (red dashed lines in Fig. 1d-e) by Eq. (6-9). Then we can get the amplitudes of zonal and meridional wind shears, which are the 5 modules of ′ ⁄ and ′ ⁄ , respectively. Finally, the GW-induced shear (S) can be calculated by Eq. (10).
A brief summary of the results from the spectral decomposition method (red) and theory (black) are shown in Fig. 1.
From Fig. 1, we can see that the GW-induced winds and shears derived from spectral decomposition method (red dashed lines) agree well with the theoretical results (black solid lines) below 100 km. The bad consistencies occur at around the upper boundary. Thus, we will focus on the results below 100 km in the following analysis. 10

GW-induced Shear Derived from SABER GW profiles and Uncertainties
A key step to derive the GW-induced shears is the extraction of GWs from the SABER temperature profile. The extraction methods of GWs from satellite data have been developed by Fetzer and Gille (1994) and improved greatly since (e.g., Preusse et al., 2002;Ern et al., 2004Ern et al., , 2011Ern et al., , 2018Chen et al., 2019;Alexander et al., 2008, 2018, Wang & Alexander, 2010Alexander, 2015). We have developed a similar method in our previous studies (Liu et al., 2014b(Liu et al., , 2017, which is summarized here. 15 First, the daily SABER temperature profiles ( ) in a latitude band of 5° are selected. Second, at each altitude, these selected data are fitted by harmonics with zonal wavenumbers ranging from 0 to 6, which are mainly planetary waves and non-migrating tides and are removed from ( ) to get the residual temperature ( ). The component of wavenumber 0 is considered as the zonal mean temperature � ( ). Third, DFT (Eq. (13)) is applied on each residual profile to amplitude ( � , ) and phase ( , ) at the vertical wavenumber ( ) with = 3, ⋯ , 18. The corresponding vertical wavelengths are from ~5 km to ~30 km. From 20 these monochromatic waves, we reconstruct a GW profile ′ ( ). The above three steps are applied on both the ascending and descending nodes, respectively, such that it minimizes the influences of migrating tides on GWs (Preusse et al., 2002).
A GW profiles pair is defined as the two adjacent SABER GW profiles, whose along-track distance is less than 400 km.
The 400 km criterion is fulfilled only for the short-distance pairs for SABER measurement ( The horizontal wavenumber ℎ , which is derived from GW profiles pair, is in the along-track direction. It is smaller than the real wavenumbers since GWs may not propagate exactly in the along-track direction. This will induce uncertainties in deriving the GW-induced shears since the GW propagation direction cannot be determined from a GW profiles pair (Ern et al., 2004). The uncertainties of the GW-induced shears are estimated in below.
The relation of the horizontal wavenumber in the along-track direction ( ℎ ) and that of in the GW propagation direction 5 ( ℎ ) can be expressed as ℎ = ℎ cos . The subscripts "a" and "r" denote the along-track and real physical quantities, respectively. The angle ( ) between the along-track direction and GW propagation direction can vary from 0 º to 360º. Here, the angle is restricted in the range of 0-90º. This restriction maps the angle of four quadrants into the first quadrant since only the magnitude of ℎ and ℎ are considered here. If ≠ 0, this will induce uncertainties of intrinsic frequencies (�) and thus GW-induced shear ( ). Here the uncertainty of S (noted by Se) is defined as ratio between the along-track (noted by ) 10 and real (noted by ). According to Eq. (10) and (4-5), we get, Equation (14) shows that if = 0 or = 0, then = 1 and the GW-induced shear is accurate. if ≠ 0 and = 2Ω sin ≠ 0, then Se > 1 and the total wind shear is overestimated. Figure 3 shows the dependencies of on horizontal ( ℎ ) and vertical wavelengths ( ), the angle and latitude . From 15 Fig. 3a, we can see that increases with the increasing ℎ and the decreasing . Figure 3(b) shows that increases with the increasing ℎ and the increasing . Figure 3(c) shows that increases with the increasing ℎ and the increasing .
Comparisons the relative importance of , ℎ , , and in changing , the angle is the most important one. If we assume that GW propagates in an arbitrary direction, is less than 1.2 (indicated by a red contour line) for a large fraction of GWs at = 30°. The fractions of GWs, which are overestimated by a factor of 0.2, increase with the increasing latitudes. This 20 indicates that the method proposed here can be used to estimate GW-induced shears even though the wave propagation direction cannot be determined from a GW profiles pair. When we analysed the derived total wind shears, the overestimation at high latitudes should be considered.

Comparisons GW-induced Shears with Model and Observational Results
Due to the uncertainties, the GW-induced shears will be further examined by comparing with models and observational results. 25 According to the wind shears during 1-10 July simulated by WACCM (Liu, 2017) and observations, we take the GW-induced shears during May-August of 2018 as an example for the purpose of comparison. The longer date coverage is chosen to cover a wider latitude range. Figure 5 shows the zonal mean and standard deviations (stds) of GW-induced shears and the top 10% largest shears during three periods of 2018. The latitude band used for zonal mean has a width of 5° with overlap of 2.5°. The latitude coverage is from 52.5°S-82.5°N or 82.5°S-52.5°N due to the yaw cycle of the SABER measurement. The three periods 30 have centers in July and extend one or two months, such that the results can illustrate the main features during summer. https://doi.org/10.5194/acp-2020-515 Preprint. Discussion started: 1 July 2020 c Author(s) 2020. CC BY 4.0 License.
Three main features can be found in the zonal mean shown in Fig. 4. Firstly, the maxima of are ~13-17 ms -1 km -1 with stds of ~9-11 ms -1 km -1 at z~90-100 km (around the mesopause) at all latitudes of the three different time intervals. Especially, the GW-induced shears reach their maxima at around 70°S. The maxima are ~10-13 ms -1 km -1 with std of ~8-10 ms -1 km -1 at z~80-90 km at latitudes higher than 40°N, where is near the summer mesopause (referred to contour lines of the zonal mean temperature). Secondly, the zonal mean shears have similar latitude-height patterns during the three different time intervals 5 although they are averaged over different time interval lengths (e.g., 31 days during 0701-0731 and 62 days during 0629-0830).

Comparisons with Model Results
The latitude-height patterns of derived here agree well those simulated by WACCM during 1-10 July (Liu, 2017).
Specifically, the GW-induced shears derived here have peaks at ~z=80-90 km at latitudes higher than 50°N during 0506-0627. 15 Moreover, there is another peak at ~z=90-100 km of the southern high latitudes during 0629-0830. These peaks agree with the WACCM simulation results that the large shears have peaks at around the mesopause (Xu et al., 2007;Fig. 2a of Liu (2017)).
However, the shear peaks derived here are at a slightly lower height than the WACCM results (Liu, 2017).
Compared to the WACCM simulation results, the differences pertain to the magnitudes of the zonal mean and specifically: (1) the maxima of the zonal mean . Figure 2(a) of Liu (2017) showed that the maxima of the zonal mean are 20 20-40 ms -1 km -1 near the mesopause at latitudes higher than 50°N/S and are ~20-25 ms -1 km -1 at latitudes lower than 50°N/S. This indicates that the GW-induced zonal mean derived here are about 70% smaller than those of the WACCM simulation results.
(2) the top 10% largest . The simulated top 10% largest (e.g., the minima of the total 10% largest ) are larger than 50 ms -1 km -1 at high latitudes and 35 ms -1 km -1 over the equator (Fig. 5 of Liu, 2017). Thus, the top 10% largest derived here (with maxima of ~30 ms -1 km -1 ) is about 60% of the simulated results. (3) the stds. The simulated stds reach their maxima of 25 ~1.4 times of the shears (Fig. 2b of Liu, 2017). The ratio of 1.4 indicates that the maxima of the simulated std is about 40 ms -1 km -1 , which is larger than those derived from observations (~12-14 ms -1 km -1 , see the middle row of Fig. 4) at latitudes lower than 50°N/S. (4) the structures of . The structures in Fig. 2b of Liu (2017) have much finer scales than those in Fig. 4. The different spatial scale coverages, which will be discussed in Sect. 5, might be responsible for the smaller values of wind shears derived from observations. 30 In general, the GW-induced shears derived here can reproduce the latitude-height pattern and 60-70% of the simulated shear magnitudes in the zonal mean sense.

Comparisons with Observational Results
To compare the GW-induced shears derived here and those observed by ground-based lidar and sounding rockets (Larsen, 2002;Yue et al., 2010), we show in Fig. 5 the profiles of and their zonal means as well as the top 10% and 1% largest during January and July at around 40°N and Equator. The January and July are representative months for winter and summer, respectively. The latitudes of Equator and 40°N may be representative for low and middle latitudes, where sounding rocket 5 measurements were performed ( Fig. 1 of Larsen, 2002;Larsen & Fesen, 2009 2010)). This agrees well with the profiles shown in Fig. 5. After averaging the shears observed by the CSU lidar during summer and winter months, the zonal mean shears are ~12-17 ms -1 km -1 with stds of ~10 ms -1 km -1 (Fig. 11c of Yue et al. (2010)). They are slight larger than the magnitudes of 10-13 ms -1 km -1 derived here (the 15 right column of Fig. 5). The zonal mean shears (the right column and upper row of Fig. 5) derived here have similar magnitudes during January and July above 90 km. However, the mean shears observed by the CSU lidar the have similar magnitudes in winter and summer at ~z=80-105 km. We note that the magnitudes of derived here are larger in July than in January at Equator. This might be related to the intra-annual variations of and will be studied in Sect. 4.
A short summary of the above comparisons is below. The GW-induced shears derived from SABER observations agree 20 with the previous observations and model results in general. This provides a global climatology of large shears around the mesopause region partially based on observations. The magnitudes of derived from SABER observations are similar to those observed by lidar and sounding rocket but are about 60-70% of the high resolution WACCM results in the zonal mean sense.
The difference probably comes from (1) the coarse horizontal samplings ~250-350 km of satellite observations and (2) only GWs with ≥ 5 km and ℎ ≥ 2∆ (∆ is the distance between the two GW profiles in a pair) used to derived shears. This 25 will be further discussed in Section 5.

Climatology of GW-induced Shears
With the advantage of 18-year SABER observations, the climatology of the GW-induced shears can be explored on the aspects of their latitudinal variations and intra-annual variations. The four seasons in the northern hemisphere (NH) are: spring (MAM: March, April, and May), summer (JJA: June, July, and August), autumn (SON: September, October, and November), winter Figure 6 shows the latitude-height contours of the zonal mean and stds of the GW-induced shears and the top 10% largest shears during four composite seasons. Each composite season is made up by the superpositions of the corresponding seasons from 2002 to 2019. The numbers of profiles used to derive GW-induced wind shears in each season (the fourth row of Fig. 6) have peaks around 50°N/S due to the changes of ascending and descending nodes. The sharp changes of the profile numbers 5 around 50°N/S might induce the discontinuity of the latitudinal variations of GW-induced shears. By further examination of the discontinuities in Fig. 6, we find that the zonal means exhibit more obvious discontinuities at around 50°N (50°S) than those at around 50°S (50°N) during spring and summer (autumn and winter). This is because there are fewer samplings (the fourth row of Fig. 6) at around 50°N (50°S) than those at around 50°S (50°N) during spring and summer (autumn and winter).

Latitudes Variations of GW-induced Shears
The hemispheric asymmetry of the sampling is induced by the inconsistency of the date coverages of yaw cycle and 10 season. Consequently, we show in Fig. 7  and summer (autumn and winter). Moreover, the wind shears have peaks at a lower height (z~80-90 km) and at latitudes of 82.5°S-50°S (50°N-82.5°N) during autumn and winter (spring and summer). These lower height peaks during spring and autumn (highlighted by blue rectangles) are weaker than those during summer and winter (highlighted by red rectangles).
Comparing with Fig. 7 (the second and third columns, 0228-0512, 0502-0715), we find that the weak peak at z~80-90 km during spring is contributed from the wind shears during May (highlighted by a red rectangle in Fig. 7 during the yaw cycle of 25 0502-0715) since there is no peak at similar location during the yaw cycle of 0228-0512. The same is true during autumn, when the weak peak at z~80-90 km is contributed from the wind shears during November (highlighted by a red rectangle in Fig. 7 during the yaw cycle 1031-114) since there is no peak at similar location during the yaw cycle of 1228-0318. The stronger peaks during summer and winter in Fig. 6 are contributed from those during the yaw cycles of 0502-0715 and 1031-0114, respectively. The stronger peaks above 90 km and at ~z=80-90 km (marked by red rectangles) are both at around the 30 mesopause as referred to the zonal mean temperature (contour lines in the second rows of Figs. 6 and 7).
The std and the top 10% largest shears, which are shown in the second and third rows of Figs. 6 and 7, respectively, have similar patterns as that of zonal mean shears. The maxima of the std and the top 10% largest shears are, respectively, ~12 ms -1 km -1 and ~30 ms -1 km -1 , which are slightly less than that shown in Fig. 4. This is because the sampling profiles in Figs. 6 and 7 (composite season or yaw cycle over 18-year) are much larger than those in Fig. 4 (only one yaw cycle in one year).
Since the patterns of zonal mean and stds of the GW-induced shears and the top 10% largest shears are similar to each other (as shown in Figs. 4, 6, 7), only the zonal mean shears during each summer from 2002 to 2019 are shown in Fig. 8. It can be seen that the latitude-height distributions of GW-induced shears, including the peaks at lower heights (around the 5 mesopause region) of high latitudes, are similar to the 18-year's mean results shown in Figs. 6 and 7. However, the GWinduced shear magnitudes (shown in Fig. 8) exhibit year-to-year variations. For example, at the SH high latitudes, the wind shears above 90 km are strongest during 2008 and 2019 and weakest during 2002. At the NH high latitudes, the GW-induced shears at ~z=85-95 km vary with year more greatly and have smaller values, as compared to those at around 80 km.

Intra-annual Oscillations of GW-induced Shears 10
Since the GW-induced shears are prominent around the mesopause region, their intra-annual oscillations will be studied at ~z=60-100 km. Figure 9 shows the monthly zonal mean GW-induced shears at four latitudes bands of the NH and SH from 2002 to 2019. A general feature of time-height variations GW-induced shears are the annual (AO) and semiannual oscillations (SAO). To quantify the exact amplitudes and phases of AO and SAO, harmonic fitting is applied on the GW-induced shears.
The fitting function has periods both AO and SAO. Figure 10 shows the amplitudes and phases of both AO and SAO at four 15 latitude bands of northern and southern hemispheres (SH).
At 50°N/S (the first row of Fig. 9 and the first column of Fig. 10), the GW-induced shears exhibit different height dependencies of AO and SAO. At 50°N, both AO and SAO reach their maxima at 80 km, while SAO has another peak at 97 km. At z=75-92 km, the AO is dominant and has peak in June. Below 75 km and above 92 km, the AO and SAO are almost equal partitioned. At 50°S, both AO and SAO reach their maxima at ~81 km, while AO has another peak at 98 km. Above 92 20 km and below 68 km, the AO is dominant and has phase in June. At z=80-90 km, the AO and SAO are almost equal partitioned and have peaks in December and June, respectively. The amplitudes of SAO at 50°N/S have similar amplitudes and height variations. However, the amplitude of AO at 50°N is smaller than that at 50°S. This makes the GW-induced shears hemispherical asymmetric.
It should be noted that the phase of AO at 50°N shifts from December at ~65 km to June at ~75 km and then shifts from 25 June at ~88 km to December at ~100 km. Whereas, the phase of AO at 50°S shifts from June at ~70 km to December at ~77 km and then shifts from December at ~85 km to June at ~65 km. In summary, in each hemisphere, the phase of AO shifts from winter below ~65 km (~70 km) to summer at ~z=75-88 km (~z=77-85 km) and then shift to winter again above 95 km at 50°N (50°S), respectively.
At 35°N/S (the second row of Fig. 9 and the second column of Fig. 10), the GW-induced shears exhibit both AO and 30 SAO. At 35°N, the amplitudes of AO and SAO vary with height in a similar pattern as those at 50°N but have smaller values.
The phases of AO and SAO are also in June when their amplitudes are prominent at z=75-92 km and then shift to winter above https://doi.org/10.5194/acp-2020-515 Preprint. Discussion started: 1 July 2020 c Author(s) 2020. CC BY 4.0 License. ~92 km. At 35°S, the amplitude and phase of AO vary with height in a similar pattern as those at 50°S. The amplitude of SAO is dominant below 90 km with peak in June.
At 20°N (the third row of Fig. 9 and the third column of Fig. 10), the amplitudes of AO and SAO exhibit similar height variations as those at 50°N and 35°N but have smaller values. The AO and SAO reach their peaks at around 81 km in June and have values of ~0.7 ms -1 km -1 and ~0.5 ms -1 km -1 , respectively. At 20°S, the SAO is in the dominant position and has peaks 5 in June at ~92 km and ~67 km. At 5°N/S and at z=85-98 km, the AO is in the dominant position and has peak shifting from July to March. At ~z=71-80, the SAO is in the dominant position and has peak shifting from March to January. of GW-induced shears agree with GWSTA and GWMF on the aspects of phase shifts and hemispheric asymmetry ( Fig. 2 and   3 of Chen et al., 2019). However, the heights at which phase shifts occur are different due to the AO and SAO in the background temperature and static stabilities (Liu et al., 2020).

Discussions
To determine the horizontal wavenumbers in the zonal and meridional directions, at least three profiles should be sampled at 20 different locations of the same wave (Wang & Alexander, 2010, Alexander, 2015, Alexander et al., 2018. For the SABER measurement, there are 15 orbits in the ascending and descending nodes, respectively. The nearest distance between two orbits is about 24°, which is much longer than the horizontal wavelengths of most of GWs. This limits our ability to deduce the zonal and meridional horizontal wavenumbers and leads to the uncertainties in deriving the GW-induced shears. The horizontal wavenumber derived from a GW profiles pair is in the along-track direction. It is in general smaller than the horizontal 25 wavenumbers of GWs in reality. The angle ( ) between the along-track direction and real GW propagation direction is the dominant source of the uncertainties. This will overestimate the GW-induced shears as shown by Eq. (14) and Fig. 3. According to Eq. (14), the influences of on increase with the increasing latitudes due to increasing inertial frequency . For the extreme case, = 1 for any since = 0 over Equator. On the other hand, at the high latitudes, GWs propagate mostly in the zonal direction since their sources are mainly jet/front, topography (Fritts & Alexander, 2003;Plougonven & Zhang, 2014). 30 Fortunately, the SABER orbit track intersects with zonal direction at a smaller angle at high latitudes than that at lower latitudes due to the changes of ascending/descending nodes. This reduces the uncertainties of the GW-induced shears for the zonally propagating waves. Thus, the uncertainties might be smaller than 1.2 as shown in Fig. 3c.
Even with an over-estimation of 1.2, the GW-induced shears derived here are smaller than those of high-resolution model simulations by a factor of 60-70% (Liu, 2017). The smaller GW-induced shears might be induced by the following two reasons: (1) the coarse horizontal samplings of satellite observations and (2) only GWs with ℎ ≥ 2∆ (Nyquist limit) and ≥ 5 km used to derived shears. The reason for (1) is that the latitude-height variations of stds in Figs. 4, 6, 7 are smoother than those in Fig. 2b of Liu (2017). The small-scale variations in Fig. 2b of Liu (2017) might be smoothed out due to the observational 5 filter of SABER limb sounding pattern (Preusse et al., 2002;Ern et al., 2018).
According to the SABER sampling ( Fig. 1 of Ern et al. 2011), the sampling distance of a GW profiles pair ∆~250 − 350 km limits the resolved GWs with ℎ ≥ 500 − 700 km. Whereas, the horizontal resolution of WACCM is about 25 km (Liu, 2017), this can resolve GWs with ℎ ≥ 50 km according to Nyquist limit, though waves with ℎ ≤ 200 km are excessively damped by numerical diffusion in WACCM (Liu et al., 2014c). The longer sampling distance eliminates the small-10 scale variations in GW-induced shears and thus reduces the magnitudes of GW-induced shears. The influence of the horizontal sampling distances on wind shears can be further confirmed by the simulation results presented by Shinagawa et al. (2017), who showed that the longitude-latitude distributions of the zonal wind shears have peak values of 16-18 ms -1 km -1 at 100 km during summer and winter (Figs. 6 and 7 of their paper). Certainly, the zonal mean of the zonal wind shears should be much smaller than 16-18 ms -1 km -1 at 100 km. This magnitude is also smaller than the simulation results presented by Liu (2017). 15 The smaller magnitude of the wind shears might result from the different resolution used by their models. The Ground-totopside model of Atmosphere and Ionosphere for Aeronomy (GAIA) used by Shinagawa et al. (2017) has a grid size of 2.8° longitude by 2.8° latitude horizontally and 0.2 scale height vertically. Whereas, the WACCAM used by Liu et al. (2014c) and Liu (2017) has a quasi-uniform horizontal resolution of∼25 km, and a 0.1 scale height vertically.
For reason (2), the cutoff criterion of ≥ 5 km is used here to get more reliable GW profile through DFT. This cutoff 20 criterion is related to the vertical resolution of the SABER measurement and is the same as that used by Ern et al. (2018) to remove the artificial oscillations. To test the influences of the cutoff criterions on the GW-induced shears, we perform a same procedure as that described in Sect. 2 and 3 but relax the cutoff criterion to ≥ 3 km. The GW-induced shears derived with cutoff criterion of ≥ 3 km are shown in Fig. 11. It can be seen that the latitude-height patterns of the GW-induced shears are the same as those shown in Fig. 4. However, the maxima of the zonal mean GW-induced shears increase from 12-17 ms -25 1 km -1 for ≥ 5 km to 15-21 ms -1 km -1 for ≥ 3 km. This illustrates that magnitudes of the GW-induced shears increase with the decreasing cutoff vertical wavelengths.

Summary
Due to the importance roles of the large shears in the MLT region on the atmospheric dynamics and electrodynamics, a method of deriving GW-induced shears is proposed in the work. The theorical basis of the method is the dispersion and polarization 30 relations of a linear GW. The data are the SABER temperature profiles measured over the past 18-year. Based on the method and the data, the global GW-induced shears are studied over a time span of 18-year.
The GW-induced shears derived here agree with previous lidar and sounding rocket observations in the aspects of height structures and magnitudes on the climatology sense. Moreover, the GW-induced shears derived here agree with the highresolution model simulation results in the aspects of latitude-height patterns but have smaller magnitudes in the zonal mean sense. The GW-induced shears reach their maxima around the mesopause region and increase with the increasing latitudes. At most latitudes and during all seasons, the maxima of GW-induced shears are at ~z=90-100 km. At high latitudes of the summer 5 hemisphere, the maxima of GW-induced shears at a lower height (~z=80-90 km). This latitude-height pattern of GW-induced shears is independent on years. The magnitudes of the GW-induced shears exhibit year-to-year variations.
The GW-induced shears exhibit more prominent AO and SAO at high latitudes than those at lower latitudes. The height variations of the amplitudes AO and SAO are hemispheric asymmetric. The strong AO occurs at around 80 km in the NH and above 92 km in the SH. At middle to high latitudes, the phases of AO shift from winter to summer and then to winter again 10 with the increasing height. The amplitudes of SAO decrease with the decreasing latitudes. The phases of SAO are in May and June when the SAO reach its peak at middle to high latitudes.
The main limitation of the method is the overestimation of the GW-induced shears due to the unresolved GW propagation direction by the method. The other limitations, such as the observational filter, long sampling distance and cutoff criterion of the vertical wavelength, will underestimate the GW-induced shears. To overcome these limitations, it is necessary to develop 15 new techniques and/or more dense satellites to observe the atmosphere with higher horizontal and vertical resolutions.
Data availability. The SABER data were downloaded from ftp://saber.gats-inc.com/Version2_0/Level2A/ Author contribution. XL and JX designed the study and wrote the manuscript. JY and HL contributed to the discussion of the 20 results and the preparation of the manuscript. All authors discussed the results and commented on the manuscript at all stages.
Competing interests. The authors declare that they have no conflict of interest. Figure 1: Synthetic temperature perturbation profiles (a, 1 ′ ( ) and 2 ′ ( ) are represented by black and red lines, respectively) and the corresponding winds (b, zonal; c, meridional) and shears (d, zonal; e, meridional). The winds and shears are, respectively, calculated from 5 the prescribed amplitudes and wavenumbers (black line, labelled as "th") and reconstructed by spectral decomposition method (red dashed line, labelled as "rc"). All panels have the same y-axis scale.