Large-scale gravity wave perturbations in the mesopause region above Northern Hemisphere midlatitudes during autumnal equinox: a joint study by the USU Na lidar and Whole Atmosphere Community Climate Model

To investigate gravity wave (GW) perturbations in the midlatitude mesopause region during boreal equinox, 433 h of continuous Na lidar full diurnal cycle temperature measurements in September between 2011 and 2015 are utilized to derive the monthly profiles of GW-induced temperature variance, T ′, and the potential energy density (PED). Operating at Utah State University (42 N, 112W), these lidar measurements reveal severe GW dissipation near 90 km, where both parameters drop to their minima (∼ 20 K2 and ∼ 50 m2 s−2, respectively). The study also shows that GWs with periods of 3–5 h dominate the midlatitude mesopause region during the summer–winter transition. To derive the precise temperature perturbations a new tide removal algorithm suitable for all ground-based observations is developed to de-trend the lidar temperature measurements and to isolate GW-induced perturbations. It removes the tidal perturbations completely and provides the most accurate GW perturbations for the ground-based observations. This algorithm is validated by comparing the true GW perturbations in the latest mesoscale-resolving Whole Atmosphere Community Climate Model (WACCM) with those derived from the WACCM local outputs by applying this newly developed tidal removal algorithm.


Introduction
Gravity wave (GW) forcing and its associated spatial and temporal spectrum within the mesosphere/lower thermosphere (MLT) are the key parameters for understanding the energy and momentum transfers between the lower and upper atmosphere GWs (Fritts and Alexander, 2003).GWs also play critical roles in various coupling processes between neutral atmosphere and the ionosphere (Vadas and Liu, 2013;Liu and Vadas, 2013).Mainly generated in the troposphere by orography, convection, and jet-front systems, the GWs propagate upward and horizontally throughout the atmosphere.Due to conservation of wave energy density, the wave amplitudes increase as they propagate into higher altitudes to compensate for the decreasing air density, until they become unstable and break mostly in MLT region.
Based on Ern et al. (2004), GW potential energy density (PED) is directly related to the GW momentum flux (MF) and the associated body forcing on the mean flow.Indeed, Yuan et al. (2016) successfully derived the momentum flux of a GW packet through this approach by using the coordinated observations from the Na lidar and a co-located mesospheric temperature mapper (Pautet et al., 2014).The calculation of PED requires precise measurements of both slow varying background temperature modulated by largescale tidal waves and the GW-induced temperature perturbations through carefully de-trending of the data.So far, experimental investigations on GW perturbations and dynamics have mostly focused on nocturnal small-scale individual events with short periods (Bossert et al., 2015;Fritts et al., 2014;Cai et al., 2014;Yuan et al., 2016).Except for a few studies in the polar region, long-period large-scale GWs, especially inertia GWs, which play critical roles in MLT dynamics, are not well studied, leaving contributions from most of these waves unchecked.The ground-based experimental studies on climatological GW perturbations (Gardner and Liu, 2007;Rauthe et al., 2006Rauthe et al., , 2008) ) have been limited to nocturnal observations as well with various de-trending algorithms developed, trying to remove tidal biases from the measurements.At the same time, several bandpass filter techniques are also be attempted and utilized in similar GW investigations.However, it is quite difficult to set the frequency range to cover the whole GW temporal spectrum, because some inertia GWs can have quite low frequencies that are close to those of tidal waves (diurnal tide excluded).For example, the GW inertial period is 18.44 h at Utah State University (USU) based on the inertial frequency ω = 2 sin φ, where is the Earth rotation rate and φ is latitude.In addition to these challenges on the reliable measurements of GW perturbations, large tidal day-to-day variability makes it extremely difficult to calculate the tidal wave modulations precisely, not to mention distilling GW perturbations, even for the spaceborne observations with global-scale coverage.
In this paper, we utilize the full diurnal cycle temperature observations (a total of 433 h) in the month of September between 2011 and 2015, taken by the USU Na lidar, to calculate the monthly temperature variance and PED during boreal autumnal equinox in the midlatitude mesopause region.
Here we develop a new algorithm to de-trend and remove the tidal bias in the lidar temperature measurements.To evaluate this new algorithm, we introduce the latest Whole Atmosphere Community Climate Model (WACCM) mesoscale simulation results (Liu et al., 2014), which generate global distribution and variations in GW.This new and unique study covers the GW spectrum up to inertia period, providing the most comprehensive large-scale GW results from groundbased observations.As the only investigation of this kind based on full diurnal cycle observations, these results could be utilized to evaluate the current effort of GWs simulations from all whole-atmosphere general circulation models, such as WACCM (Garcia et al., 2007), that can simulate GW impacts and evolution from the surface up to the lower thermosphere.

Instrumentation and model involved
The USU Na lidar is a narrowband resonance fluorescence Doppler lidar system operating at the Na D 2 line with a 120 MHz full width at half maximum (FWHM) laser pulse bandwidth.It has the capability of measuring highresolution temperature and horizontal wind profiles within the mesopause region (80-105 km) in full diurnal cycles (Krueger et al., 2015).Collaborating with the other surrounding instruments, such as the Advanced Mesospheric Temper-ature Mapper (Pautet et al., 2014) and Meteor Wind Radar, the lidar can provide critical measurements of tidal and planetary wave variations (She et al., 2004;Yuan et al., 2008Yuan et al., , 2010)).The high-resolution lidar measurements also provide detailed description on GW activities and their impacts in the MLT (Yuan et al., 2014(Yuan et al., , 2016;;Cai et al., 2014;Fritts et al., 2014).
WACCM is a comprehensive numerical model covering the altitude range between the surface and thermosphere.It uses the National Center for Atmosphere Research (NCAR) Community Earth System Model (CESM) as a common numerical framework (Garcia et al., 2007).Liu et al. (2014 and the references within) developed a high-resolution WACCM run based on the framework of WACCM version 5 and simulated GW activity in mesoscale globally from the troposphere to the lower thermosphere.
The lidar data are first processed with 1 h temporal resolution and smoothed vertically through a moving 2 km FWHM Hanning window.The temperature measurements with uncertainty equal or larger than 15 K are treated as bad data.There are a total of 433 h of continuous full diurnal cycle lidar data at USU involved in this study.These September data include 116 h in 2011; 73 h in 2012, with one gap of less than a day; 64 h in 2013; 89 h in 2014; and 91 h in 2015.The WACCM temperature simulations included in this study are from 1 to 10 September at 41 • N from 0 to 360 • longitude, with zonal resolution of 0.3125 • and temporal resolution of 1 h.As with the lidar data processing, the WACCM hourly temperature outputs are also smoothed with a 2 km FWHM Hanning window in the vertical direction.

Lidar data
From Ern et al. (2004), the PED can be calculated as follows: where g = 9.5 m s −2 is the gravitational constant; T is the mean temperature with tidal modulations removed; T 2 is the variance of the temperature perturbation T , which is the fitting residual that will be discussed later; and is the background Brunt-Vaisala frequency squared with the lapse rate g C P = 9.5 K km −1 in mesopause region.
To calculate GW-induced temperature perturbations at each lidar-measured altitude, the slowly varying background temperature, T 0 , has to be calculated and subtracted from the lidar data.This can be achieved through a two-step algorithm.First we apply a least-squares fitting (LSF) algorithm (shown in Eq. 2) that covers all the tidal periods on the multi-day continuous lidar measurements to derive the mean temperature and tidal period perturbations, provided that tidal period inertia GWs do not exist or are insignificant.Second, by utilizing the derived mean temperature and tidal period perturbations in step 1 we reconstruct the background temperature, T 0 , in the mesopause region as follows: where T is the mean temperature, i = 1 for the diurnal tide, i = 2 for the semi-diurnal tide, i = 3 for the terdiurnal tide, i = 4 for the quardiurnal tide, and φ i represents the tidal phase for the ith tidal component.Considering the large dayto-day tidal variability, a 24 h fitting window sliding across the multi-day lidar measurements at 1 h step is applied.The abovementioned LSF is conducted within each of the fitting windows to derive the tidal period perturbations that are regarded as the representatives of solar thermal tidal waves within the fitting window.These derived tidal period perturbations, along with the simultaneously deduced mean temperature, T , are utilized to reconstruct the slow varying background temperature, T 0 .After de-trending the temperature at each lidar altitude, the differences between the lidar temperature measurements and this reconstructed background temperature are treated as GW-induced temperature perturbations, T .The temperature variance, T 2 , for PED calculation shown in Eq. ( 1) can be calculated within the fitting window during the same process.The window then moves 1 h forward to repeat the above process until it reaches the end of the data set.This algorithm is applied on the full diurnal cycle lidar temperature measurements for each September to achieve the monthly profiles of temperature variance and PED of the year.The final mean monthly profiles of these two parameters, representing the mean GW activities, are achieved by averaging all these September results.Figure 1 shows an example of this new method of determining the GW-induced temperature perturbations based on a Na lidar campaign from UT day of year (DOY) 252 to 255 (9 to 12 September) in 2015.Figure 1a is the comparison between the hourly lidar temperature measurements and the reconstructed slow varying background temperature, T 0 , at 89, 92 and 95 km based on the LSF described in the previous paragraph.Figure 1b is the corresponding T at these three altitudes, which is the difference between the lidarmeasured temperatures and the reconstructed T 0 .As shown in the Fig. 1a, the reconstructed T 0 from this new method fits the lidar temperature measurements extremely well, reproducing the large-scale day-to-day variability in the MLT temperature.It also proves that the new algorithm generates realistic background temperature and GW-induced temperature perturbations.
Figure 2 indicates the monthly averaged temperature variance profiles (Fig. 2a) and PED profiles (Fig. 2b) derived from this LSF algorithm.On the other hand, to investigate the differences between these full diurnal cycle results and those based solely upon nocturnal lidar measurements, we reprocess the USU Na lidar nighttime data between 02:00 and 12:00 UT by following the algorithm presented by Gardner and Liu (2007).In this algorithm, the lidar-measured temperature variations at each lidar altitude during the whole night are linearly fitted to a straight line to calculate the background temperatures, T 0 .The justification is that large-scale waves, such as tides and planetary waves, are all slowly varying and can be treated as constant during the night.T 0 is later subtracted from the lidar-measured temperature to retrieve temperature perturbations.This is followed by another linear fit in the vertical direction on the resulting temperature perturbation profiles to calculate and remove the potential tidal biases residual in the first fitting attempt.The justification is that the tidal waves, especially the dominant semidiurnal tide at midlatitudes, usually have very long vertical wavelength compared to the lidar range.Thus, their effects could be further removed by this second linear fitting.These resulting profiles from this approach, termed nightly linear fit (NLF), are presented in Fig. 2 as well, alongside with those derived from the LSF on the lidar's full diurnal cycle data.

WACCM data
To evaluate this new LSF algorithm described above and see whether the deduced GW perturbations can precisely represent those of the whole large-scale GW temporal spectrum, we introduce the latest high-resolution WACCM simulations in September that are able to simulate global GW distribution and variations.Here the WACCM data are processed in two completely different approaches and two sets of GW perturbation results are generated.The first approach is to process WACCM's high-resolution outputs at the USU location with the same LSF tidal fitting algorithm discussed above, generating background T 0 and GW-induced temperature perturbations T and variances T 2 along with the associated PED profiles.Again, we also conduct the two-step fitting algorithm (NLF) for nighttime WACCM data at USU to study the potential differences.The second approach is to fit each hour's WACCM temperature data from longitude 0 to 360 • at each altitude at 41 • N with global-scale wave components of zonal wavenumbers 0-10.These removed components with zonal wavenumber 0-10 are most likely induced by the major large-scale waves, such as tidal waves and planetary waves, including the large-scale waves due to tidal-planetary wave interactions.The background temperature T 0 can be reconstructed with these low wavenumber components, along with the zonal mean temperature.The fitting residual of this zonal wavenumber removal approach, named as ZW10R for the rest of the paper, are treated as GW-induced temperature perturbations, T .Finally, the T 0 , T (zonal mean temperature) and T at USU location are chosen to calculate T 2 and PED profiles, as shown in Fig. 3b and c.This algorithm distinguishes GWs and tidal/planetary waves from their spatial scale.Thus, with the complete 24 h coverage in this WACCM simulation, the associated results should be the most reliable representation for the GW-induced temperature perturbations and the related body forcing results generated within this model simulation.

Results and discussion
As shown in Fig. 2a and b, the derived lidar temperature variance and PED profiles for September share very similar vertical structure within the mesopause region between 84 and 99 km.The errors represent the averaged goodness of fit (chi-square divided by the difference between the numbers of sampling and fitting parameters) in the temperature variance, T 2 profiles and those propagating into the PED calculation.Both slowly decrease to their minima from 84 km up to ∼ 89-90 km and then increase quickly above this altitude.For example, the GW-induced temperature variance drops from ∼ 50 K 2 at 84 km to ∼ 20 K 2 near 90 km but starts to increase fairly quickly going into higher altitudes, close to 90 K 2 near 99 km.It is worth mentioning that such a "node" feature in temperature perturbation within mesopause region has been reported in nocturnal lidar observations at another midlatitude location (Rauthe et al., 2006) and lowlatitude station (Lu et al., 2009).This lidar-observed "node" feature suggests a strong GW dissipation/saturation region near/below 90 km, where GWs deposit their momentum and energy to the mean flow.This statement is partially supported by large turbulence between 65 and 90 km observed by the Middle and Upper atmosphere Radar (MU Radar) (Tsuda, 2014).This feature could be interpreted by the simplified relation between the GW horizontal wind perturbations, u , and its vertical wavenumber, m, proposed by Booker and Bretherton (1967), which indicates dramatic increasing of u and m when the GW is approaching its critical level.However, to our knowledge no theoretical simulation has been able to duplicate these observations, and therefore no solid dynamic scheme has been provided.
Although T 2 and PED from WACCM do not have minima at 90 km as in the observations, they do show local minima in the upper mesosphere (∼ 82 km in ZWR10 and ∼ 78 km in LSF, respectively).Based on GW linear saturation theory, when approaching its critical level, the GW becomes unstable and experiences saturation, leading to damping of GW amplitude due to the dissipation process (Fritts and Alexander, 2003).Therefore, the "node" or the local minimum feature is likely associated with the GW critical level decided by the local horizontal wind profile.Figure 3d shows the comparison between the lidar-measured mean zonal wind in September and the zonal mean zonal wind profile on 5 September in WACCM.As shown in the figure, although the two wind profiles are both eastward in the mesopause region, they are quite distinct.The lidar-measured zonal wind is mostly slowing down as the altitude increases, while the zonal mean zonal wind in WACCM reverses from westward to eastward near 79 km and is accelerating eastward below 90 km.It is worth noting that, in the mesopause region below 88 km, the lidar mean zonal wind is much larger than the zonal wind of the model.The difference in background wind and wind reversal level, which is most likely due to deficient GW forcing in this simulation, probably results in the altitude difference of GW filtering between the simulation and observations.This difference prevents direct comparisons between the model results and those based on the lidar measurements.
Compared to the LSF results with those based on nocturnal data using NLF algorithm, there are considerable differences below around 90 km in the lidar results, shown in Fig. 2. For example, near 86 km, the results deduced from NLF are ∼ 50 % more than those derived from LSF.This is most likely due to the bias from the diurnal tide, which peaks during equinox at midlatitudes (Yuan et al., 2010), that could not be completely removed in the NLF algorithm.The figure also shows that, although NLF results are still mostly larger than LSF results, the differences between the two are much less above 90 km.When we apply the two algorithms on WACCM local outputs, the two sets of WACCM profiles (Fig. 3b and c) show even larger differences: near 75 km and between 80 and 92 km, the results from NLF overestimate the GW perturbations by almost or more than 100 %.The significant differences between 80 and 92 km in WACCM results are, again, most likely related to the bias caused by diurnal tide, while those near 75 km could be induced by planetary waves.Indeed, it has been reported that transient planetary waves could be quite active during autumnal equinox (Liu et al., 2007).Although not completely removed in the LSF tidal remove algorithm, the planetary wave effects are decreased considerably in the 24 h fitting window.
To validate the new LSF algorithm and the associated GW results above, we present the two kinds of profiles derived from WACCM data using the two completely different algorithms mentioned in the previous Section.Figure 3a shows the mean temperature T calculated from LSF algorithm at single location is mostly the same as the mean temperature derived globally from the ZWR10 approach.As indicated in Fig. 3b, the two T 2 profiles are also very similar and the differences between the two are mostly within the fitting uncertainties.Within the mesopause region, both profiles indicate the temperature variance is less 20 K 2 near 80 km and grow quickly and continuously above, as large as ∼ 50 K 2 near 95 km.However, above ∼ 97 km, the two become quite different: those deduced from ZW10R are considerably larger (close to 70 K 2 at 100 km) than those calculated from local LSF at the same altitude (∼ 40 K 2 at 100 km).One possibility is the presence of GWs with periods near 6 and 8 h, whose amplitudes grow significantly above this altitude.These waves would be excluded from the LSF timefitting method but will remain in the ZWR10 spatial filtering method.The associated PED profiles in Fig. 3c behave similarly to their related temperature variance profiles except that, above ∼ 97 km, the differences LSF and ZWR10 results are all within the fitting uncertainties in the two PED profiles.Intriguingly, there is no "node" structure in both the LSF profiles and ZW10R profiles.Overall, based on the comparison between the WACCM GW results derived from LSF and ZW10R, both algorithms generate very similar GW temperature perturbations and the PED.Therefore, both can be utilized to investigate the GWs with high confidence.For the ground-based studies without horizontal coverage, the LSF tidal removal algorithm on full diurnal cycle measurements provides the most reliable GW results.
Figure 4 illustrates the Lomb spectrum power for the temperature perturbations during each of the lidar campaigns in September between 2011 and 2015.The lidar temperature perturbations are derived from the aforementioned LSF tidal removal algorithm.Here, only modulations with significance level larger than 50 % are shown.The figure shows the dominance of long-period large-scale GWs' modulations in the mesopause region with periods between 3 and 5 h during almost all of the six campaigns, except the campaign during DOY 251 and 253 in 2012.This indicates that the GW perturbations in the mesopause region are mostly generated by this part of the GW spectrum during autumnal equinox around midlatitudes.To confirm this conclusion, we re-analyzed the lidar data with 15 min resolution and conducted the same LSF algorithm.This extends the highest GW frequency in this study from the current 0.5 to 2 h −1 , covering more highfrequency components in the total GW spectrum.The Lomb-Scargle results based on these high-resolution data again demonstrate the consistent dominance by this part of GW spectrum in midlatitude mesopause region during this part of the year.

Summary and conclusion
By applying a newly developed LSF tidal removal algorithm to the unique full diurnal cycle Na lidar observations at USU in September between 2011 and 2015, the monthly averaged profiles of temperature variance and the associated PED induced by the large-scale GWs are derived in the mesopause region during the boreal autumnal equinox near midlatitudes.The study covers the GW spectrum from 2 h up to the inertia period, providing the most comprehensive large-scale GWs results in MLT.It reveals a "node" structure near the middle of mesopause region in both profiles, decreasing GW modulations between 84 and 90 km, but a reversed trend above 90 km, where temperature variance increases from its minimum of less than 20 K 2 up to over 90 K 2 near 99 km.This "node" feature indicates a GW dissipation layer near the middle of the mesopause region that cannot be resolved or well presented in current general circulation models (GCMs).The possible mechanism may be related to the GW critical level near 90 km that prevents some GWs propagating beyond this altitude, while the GWs at other parts of the spectrum could penetrate or leak through the layer due to their fast horizontal phase speed.The dominance of long-period GWs in MLT region, especially those with periods between 3 and 5 h, is evident in the lidar observations, indicating that these GWs likely contribute the most to the GW energy and momentum in the MLT region.On the other hand, the results based solely on nighttime observations and NLF approach are found to overestimate the GW perturbations by ∼ 50 % near and below 90 km due to the bias of diurnal tide that cannot be removed completely in this algorithm.A similar test on the WACCM local results indicates an overestimate the GW perturbations of close or more than 100 % near the 75 and 80-90 km range, using NLF on nocturnal data alone.Thus, full diurnal measurements are critical for precise calculation on the GW perturbations.
To validate this new algorithm for de-trending lidar measurements, the latest high-resolution simulation results from NCAR WACCM are introduced.Besides applying the same LSF tidal removal approach on the WACCM data at the USU location, we also apply a spatial zonal wavenumber removal algorithm to the WACCM data of all longitudes at 41 • N to remove all tidal/planetary-scale modulations with zonal wavenumber less than and equal to 10.The fitting residuals from the two completely different algorithms are treated as the GW-induced temperature perturbations.The investigation reveals that, within fitting uncertainties, these two sets of WACCM GW temperature variance and PED profiles are almost identical in the MLT.This comparison clearly demonstrates that the new LSF tidal removal approach generates the most accurate GW results for the ground-based observations.This study builds a concrete foundation for future investigation on the climatology of GW perturbations and forcing for other ground-based observations.

Data availability
The lidar data of this study are available at the CRRL Madrigal database (2017) at http://madrigal.physics.colostate.edu/htdocs/.

Figure 1 .
Figure 1.(a) The lidar temperature measurements (black) and the reconstructed background temperature (red) at 89, 92, and 95 km from UT day 252 to 256 in 2015.(b) The corresponding temperature perturbation.The temperatures at 92 and 95 km are increased by 50 and 100 K, respectively, and perturbations at 92 and 95 km are increased by 20 and 40 K, respectively.The magenta lines represent the zero-point line for the temperature perturbation.

Figure 2 .
Figure 2. The 5-year USU Na lidar September monthly climatology of (a) temperature variance and (b) potential energy density (PED) obtained from the LSF method (black) and the linear fit for nightly linear fit (NLF) (red).

Figure 4 .
Figure 4. Lomb spectra power of the lidar-measured temperature perturbations during each of the six lidar campaigns in September between 2011 and 2015.The DOY of each campaign is listed in the parentheses.Only components with significant level larger than 50 % are plotted.