Seasonal Variation of Vertical Heat and Energy Fluxes due to Dissipating Gravity Waves in the Mesopause Region Over the Andes

Vertical energy transports due to dissipating gravity waves in the mesopause region (85–100 km) are analyzed using over 400 h of observational data obtained from a narrow‐band sodium wind‐temperature lidar located at Andes Lidar Observatory (ALO), Cerro Pachón (30.25°S, 70.73°W), Chile. Sensible heat flux is directly estimated using measured temperature and vertical wind; energy flux is estimated from the vertical wavenumber and frequency spectra of temperature perturbations; and enthalpy flux is derived based on its relationship with sensible heat and energy fluxes. Sensible heat flux is mostly downward throughout the region. Enthalpy flux exhibits an annual oscillation with maximum downward transport in July above 90 km. The dominant feature of energy flux is the exponential decrease from 10−2 to 10−4 W m−2 with the altitude increases from 85 to 100 km and is larger during austral winter. The annual mean thermal diffusivity inferred from enthalpy flux decreases from 303 m2 s−1 at 85 km to minimum 221 m2 s−1 at 90 km then increases to 350 m2 s−1 at 99 km. Results also show that shorter period gravity waves tend to dissipate at higher altitudes and generate more heat transport. The averaged vertical group velocities for high, medium, and low frequency waves are 4.15 m s−1, 1.15 m s−1, and 0.70 m s−1, respectively. Gravity wave heat transport brings significant cooling in the mesopause region at an average cooling rate of 6.7 ± 1.1 K per day.

where κ = R/C p , R is the dry air constant, c p is the specific heat at constant pressure, p is the atmospheric pressure, and p 0 is a reference pressure typically defined as 1,000 hPa, 〈w′θ′〉 can be obtained by linearizing θ as: 〈w′p′〉 is defined as the vertical kinetic energy flux induced by GWs. In the MLT, GWs mostly propagate upward so the overall 〈w′p′〉 is positive. As mentioned earlier, 〈w′θ′〉 = 0 for freely-propagating GWs therefore 〈w′T′〉 > 0 in the MLT. For strongly dissipated GWs, both enthalpy and sensible heat fluxes become negative. They are identical when the Boussinesq approximation is applicable (p′ is negligible and 〈w′p′〉 ≈ 0), which is generally not the case in the mesopause region. Since a Na lidar measures T′, not θ′, only the sensible heat flux 〈w′T′〉 can be directly calculated. However, if the energy flux 〈w′p′〉 can be estimated, it is possible to determine 〈w′θ′〉 from Equation 2.
Previous studies of the GW fluxes with lidar used observations at SOR and Maui (Gardner & Liu, 2007;A. Z. Liu & Gardner, 2005). The main purpose of this paper is to quantify GW fluxes 〈w′T′〉, 〈w′p′〉 and 〈w′θ′〉 in the mesopause region over the midlatitude Andes based on Na lidar data acquired in the past few years. The Andes is a hight mountain region and is different from SOR and Maui, thus has very distinguished wave generation mechanisms and background atmosphere. This work is therefore a valuable addition to the two sensible heat flux measurements at Maui and SOR.
More importantly, this is the first time that GW enthalpy flux 〈w′θ′〉 and kinetic energy flux 〈w′p′〉 are inferred from a Na lidar measurement. The results enable us to better quantify the GW dissipation since 〈w′θ′〉 is exactly zero for non-dissipation GWs while 〈w′T′〉 is not. The kinetic energy flux 〈w′p′〉 also provides additional insight to the GW propagation and energy transport in this region. As mentioned earlier, GUO AND LIU 10.1029/2020JD033825 2 of 16 quantifying the vertical fluxes is challenging because they are small statistical parameters and subject to large uncertainties. Therefore, error analysis is crucial for reliable estimate and accurate interpretation. Multiple sources that contribute to the measurement errors, including photon noise and statistical uncertainty of 〈w′T′〉, have been analyzed in the previous study by GL07 and will not be repeated here. In this study, we will analyze two additional sources of error associated with the sampling frequency and integration time using a numerical method to provide a more complete description of the heat flux error. The analysis is presented in Appendix A.
The paper is organized as follows. Section 2 describes the Na lidar data set and the method applied to derive GW perturbations. Section 3 presents the algorithm for calculating energy flux from GW temperature spectra. In Section 4, we present the seasonal variations of heat, enthalpy, and energy fluxes, the derived thermal diffusivity, and the associated heating/cooling rates. In Section 5, the contributions of GWs at different temporal scales to these fluxes are presented. The conclusion is given in Section 6.

Na Lidar Data
The observational data set was acquired with a narrow-band resonance-fluorescence Na lidar located at Andes Lidar Observatory (ALO) (30.25°S, 70.74°W) in Cerro Pachón, Chile. The lidar performs nightly measurements of temperature, line-of-sight (LOS) winds, and Na density profiles between 80 and 105 km. The ALO lidar system is equipped with four 75-cm diameter telescopes pointing toward zenith (Z) and 20° off-zenith to the north (N), south (S), and east (E) directions. In May 2014, the system was upgraded by replacing the old Coherent Ring Dye Laser with a high-power amplified diode laser (TA-SHG from Toptica Photonics) as the master oscillator (A. Z.  and the receiver system was improved with a more efficient optical design . These improvements result in a much more reliable and stable system with several times higher signals, which makes it possible to achieve many more nights of measurements with better data quality and higher resolutions. The ALO lidar has been operated regularly in terms of campaigns each year since the upgrades. This study utilizes the available data between 2014 and 2017, and the data cover every calendar month for seasonal analysis. During different campaigns, several laser beam pointing modes were used, including zenith only, zenith and off-zenith. The modes were named according to the telescopes pointed to different directions. For example, the mode ZSE means measurements were acquired in a sequence of zenith (Z), south (S), and east (E). In this study, our focus is on the vertical fluxes so only zenith measurements were selected, which is about 1/2 (in ZSZE mode) or 1/3 (in ZSE mode) of the total operation time. In order to perform seasonal analysis, we calculate the fluxes for each calendar month, using data not only in each month, but also data during the half month before and after, to increase the amount of data and reduce uncertainties. For example, the fluxes calculated for January included data from 16 December to 14 February. The number of observational nights and zenith data hours obtained in each month are listed as rows one and two in Note. The data used in each month include observations in the previous and the next half months.

Table 1 Numbers of Nights and Hours of Zenith Lidar Measurements Acquired Between May 2014 to December 2017 at ALO (Rows 1 and 2) and Total Nights and Hours Used for Heat Flux Calculation (Rows 3 and 4) in Each Calendar month
For comparison, among the previous heat flux studies with lidar, the total hours of lidar data at SOR was 370 in GL07 and 250 at Maui in A. Z. Liu and Gardner (2005), with about one-third being zenith measurements. Therefore, the ALO zenith data of 459.4 h is several times higher than the SOR and Maui data. However, only the annual mean profile was obtained from Maui data and only the semi-annual and annual harmonics of seasonal variation were estimated from SOR data. The ALO lidar data presented here is capable of resolving variations at time scales as short as two months.
To derive the LOS winds and temperature, raw lidar photon count data is normally processed at 60-s (90-s on some nights) integration time and 500-m range resolution. At these resolutions, the typical root mean square (rms) errors for temperature and LOS winds due to photon noise are respectively 1.4 K and 1.1 m s −1 at the peak of Na layer and 2.2 K and 2.0 m s −1 at 85 and 100 km. These errors increase quickly beyond this altitude range. Thus our flux calculations are limited in the height range between 85 and 100 km, where Na density is high and measurement errors caused by photon noise are relatively small. To calculate the vertical fluxes, we re-processed the data to derive temperature and wind with an increased range resolution of 2 km (but still at 500-m interval by using a binning method) to further reduce the photon noise error by about 50%. This reduction in vertical resolution filtered out GWs with vertical wavelength less than 2 km, but has negligible effect on the total heat flux because the heat flux is mostly contributed by GWs with large vertical wavelength (hence large w′). Reducing the vertical resolution while maintaining the 1-min temporal resolution is an optimal choice for including most GWs contributing to heat flux while minimizing photon noise errors.

GW Perturbations
To extract GW perturbations at well-defined frequency ranges, we used low-pass filters at four cut-off frequencies to successively separate the data into five frequency bands, with periods (1) longer than 6 h, (2) 3-6 h, (3) 1-3 h, (4) 5 min to 1 h, and (5) < 5 min. For example, frequency band (2) was obtained by subtracting 6 h low-pass filtered perturbations from 3-h low-pass filtered perturbations. We consider (1) as the background determined by long time-scale variations such as tides, and (5) as oscillations at frequencies higher than GWs.
(2), (3) and (4) are considered as low-frequency, medium-frequency, and high-frequency GWs, respectively. Although there is no clear cut in frequency separation between tides and GWs, and the 5 min is not a definitive cut-off between GW and non-GW perturbations (such as acoustic waves) due to Doppler effect of the background wind, this is a reasonable choice that includes most GWs detectable by the lidar and almost all fluxed induced by GWs. Separation of low-, medium-, and high-frequency GWs are also subject to overlaps when applying low-pass filters to cutoff frequencies. The sensible heat flux was then calculated for the total GW perturbations and individually for the 3 GW frequency bands.
As described in Section 2.1, the perturbations are grouped according to calendar months, with each month including data that span a two-month period. For each group, the fluxes (covariances) and variances are calculated at each altitude (500 m interval), using data that span 1.5 km (including data 500 m above and below). This means the fluxes and variances are calculated at an effective 1.5-km vertical resolution (but still 500-m interval). This approach further reduces uncertainties, but does not cause any loss of information because the original temperature and wind resolutions are 2 km.

Energy Flux
Gravity wave energy flux 〈w′p′〉 has rarely been examined in the MLT because of the challenge in measuring pressure and vertical velocity simultaneously and continuously over tens of hours. Only a few numerical model simulations (e.g., Ehard et al., 2016;Hickey & Cole, 1988;Yu & Hickey, 2007), and some limited observations (e.g., Gossard, 1962;A. Z. Liu, 2009;Vincent, 1984) contributed to our understanding of GW induced energy flux. It is acknowledged that this energy flux cannot be neglected and must be taken into consideration when evaluating the effects of GWs (Becker & Schmitz, 2002;Becker, 2004;Gardner, 2018;. To estimate the GW energy flux, A. Z. Liu (2009) derived an approximated expression of 〈w′p′〉 based on GW polarization relations using a broad wave perturbation spectrum that lidar can observe. This approximation assumes that all GWs propagate upward and contribute to positive upward energy flux. Although in GUO AND LIU 10.1029/2020JD033825 4 of 16 the mesopause region gravity waves are primarily upward propagating, downward propagation also exists when waves experience reflection or being ducted. Thus energy flux is likely overestimated when downward propagating GWs are ignored. Recently Gardner (2018) proposed a modification to take into account both upward and downward propagating GWs and employed a more realistic non-separable spectrum for the GW temperature fluctuations (see his Equation 32). If fraction of downward propagating wave energy is α dwn then the upward energy flux is reduced by a factor of 2α dwn . α dwn was assumed to be 0.15 based on a previous study by Hu et al. (2002) that showed that 84.4% of GWs were upward propagating in the mesopause region. Here we adopt this approach by Gardner (2018) to calculate energy flux with the ALO lidar observations.
First the frequency (ω) and vertical wavenumber (m) spectrum of gravity waves are computed with the nightly temperature perturbations. Slopes are obtained through a linear regression fitting under the logarithmic coordinates after the noise floor was removed. Shown in Figure 1 are the monthly slopes of ω and m power spectra at ALO. The annual mean slopes are −1.90 ± 0.14 and − 2.75 ± 0.31 for ω and m spectra respectively. Note that the ω spectrum has smaller uncertainties than the m spectrum because of more available data points in the temporal dimension. Meanwhile the small altitude range (15 km) and low vertical resolution (2 km) limits the accuracy in fitting the m spectrum as well as the range of vertical scales lidar data can resolve. Additionally, the m spectrum may vary with larger-scale influences and background fields. The slopes are consistent with the results reported in Gardner (1998) for SOR, which are − 1.98 for ω spectrum and − 2.74 for m spectrum, as well as other observations in the mesopause region (e.g., Lu et al., 2015b;Senft & Gardner, 1991;Yang et al., 2008).
However, the method to estimate energy flux by A. Z. Liu (2009) and Gardner (2018) needs the slope of intrinsic ω spectrum while the ω spectrum described above is non-intrinsic, as the case for all groundbased instruments. Doppler shift by horizontal winds results in differences between intrinsic and observed frequencies and as much as ±10%-20% difference in the spectral slopes depending on whether GWs are predominantly propagating along or against the background winds. This will lead to under or overestimated energy flux. For the range of wave period (5 min-6 h) and vertical wavelength (2-15 km) considered in this study, this bias can be estimated according to Equation 32 in Gardner (2018). A 15% difference between the slopes of observed and intrinsic ω spectra corresponds to 24% difference in energy flux. We regard this as the upper limit of estimated energy flux bias, because the GWs are generally propagating in multiple directions as often demonstrated airglow imager measurements (e.g., Li et al., 2011, for Maui). The actual bias is expected to be much smaller.
Presented in Figure 2 are the monthly temperature variance and vertical energy flux in the mesopause region. The temperature variance is a measure of GW intensity. There is a strong annual variation of the GUO AND LIU 10.1029/2020JD033825 5 of 16  temperature variance where the maximum ∼ 90 K 2 occurs in July and the minimum ∼40 K 2 is in December. This shows that GW are stronger during austral winter solstice and weaker during austral summer solstice. The dominant feature of energy flux (Figure 2b) is the exponentially decrease over two orders of magnitude from 10 −2 Wm −2 to 10 −4 Wm −2 with increasing altitude. Seasonally, the maximum and minimum energy fluxes are observed in August and March, respectively.

Seasonal Sensible Heat and Enthalpy Fluxes
As mentioned earlier, characterizing vertical sensible heat flux is difficult because it is a small quantity and natural variations of instantaneous flux are large. Therefore, the uncertainty of measured sensible heat flux should be addressed carefully. Assuming the measurement errors caused by photon noise for T and w are ΔT and Δw, we can express the measured heat flux as where the second term on the right hand side is the bias caused by photon noise. The upgrade of ALO lidar in May 2014 resulted in several times higher signal levels. At 60 s and 2 km resolutions, the raw photon count at peak frequency can exceed 20,000 and the corresponding bias is less than 0.02 Km s −1 , which is negligible compared with the typical magnitude of 〈w′T′〉 at 1-3 Km s −1 . For lidar systems with much lower signal levels, this bias can reach a few 0.1 Km s −1 with relative error as large as 50% and affecting the heat flux results substantially.
Apart from the bias caused by the correlation of lidar measurement errors, the nonlinearity of photomultiplier tubes (PMT) used in all lidar systems would also introduce bias to the vertical flux calculation. This is especially significant for systems with high signals. This bias can be removed if the nonlinearity of the PMT is known, but laboratory measurement of PMT nonlinearity is not accurate enough for this purpose. At ALO, this bias has been removed following a novel calibration processed based on measured temperatures at different signal levels with the same PMT as described in A. Z. . Without this careful calibration sensible heat flux calculation is susceptible to very large bias.
Another potential error source is the sampling error resulting from the limits of data resolutions, especially when high-frequency GWs are considered as in this study. Because the sensible heat flux is a result of small departure of the phase difference between w′ and T′ from π/2, ambiguity in resolving phases of the perturbations can potentially affect the sensible heat flux calculation. A detailed analysis of the sampling error is presented in the appendix. The analysis shows that this error causes a reduction of sensible heat flux of about 12% for 5-min waves and 1% for 20-min waves. As expected, this error mainly affects very high-frequency GWs close to the buoyancy frequency and decreases very quickly as wave period increases. Even though higher frequency GWs tend to contribute more to the sensible flux, this error is still very small compared to the total heat flux.
The seasonal and altitudinal variation of the sensible heat flux 〈w′T′〉 is shown in Figure 3a. The enthalpy flux is shown in Figure 3b. Note that the magnitude of 〈w′θ′〉 depends on the magnitude of θ which is in turn dependent on the reference altitude or the reference pressure p 0 in Equation 1. If we were to use the standard definition p 0 = 1000 hPa, θ and 〈w′θ′〉 would be very large. To make it easier to compare with the sensible heat flux, we used p 0 = 0.19 Pa that corresponds to the typical pressure at 90 km. Therefore, only the signs and variations of magnitude are relevant in Figure 3(b). It is obviously that the enthalpy flux is downward throughout the region and year. Strong transport happens in several months mostly in the upper region. This downward flux has a clear seasonal variation that peaks in July. A second downward flux peak is in February. The downward flux also increases with altitude, indicating stronger GW dissipation at higher altitudes. The results are remarkably consistent with our understanding that dissipating GWs always induce a downward 〈w′θ′〉 on average. The stronger dissipation at higher altitude can be expected as more waves dissipate when their amplitudes become larger. The peak dissipation in July is also consistent with the notion that strong mountain waves are present mostly in austral winter and experience strong dissipation in the mesopause region.
The wave heat transport can be described by the thermal diffusivity k H inferred from vertical enthalpy flux as Shown in Figure 4 are the monthly and annual values of k H as a function of altitude. Same as the enthalpy flux, k H is largest (≃ 800 m 2 s −1 ) in July when GW dissipation is strongest (Figure 4a). Small k H in April below 88 km and December are mainly due to weak GW dissipation. The annual thermal diffusivity inferred from enthalpy flux decreases from 303 m 2 s −1 at 85 km to minimum 221 m 2 s −1 at 90 km then increases to 350 m 2 s −1 at 99 km ( Figure 4b). Both the peak and mean values of k H at ALO are comparable to that obtained at SOR (A. Z. Liu, 2009), confirming this large heat transport by GWs.
The net thermal effects of dissipating GWs include the heating from irreversible conversion of mechanical wave energy (e.g., Gavrilov, 1990;Hines, 1965;Medvedev & Klaassen, 2003;Pitteway & Hines, 1963) and a cooling by the downward transport of sensible heat (e.g., Gardner & Liu, 2007;Hickey et al., 2011;Schoeberl et al., 1983;Walterscheid, 1981). Not all mechanical wave energy is converted to heat as it could be GUO AND LIU 10.1029/2020JD033825 7 of 16  converted into the kinetic energy of the mean winds as well (Hickey & Brown, 2000). Therefore, the heating associated with mechanical wave energy process is complicated and will not be discussed here. The dynamical cooling effect due to downward heat transport can be characterized as the heating rate calculated from vertical gradient of sensible heat flux (HF) as where ρ is the atmosphere density.
The seasonal variation and annual mean of this dynamical heating rate are presented in Figure 5. Monthly heating rate is calculated from the seasonal sensible heat flux 〈w′T′〉 (Figure 3a). The strong cooling occurs around austral winter between 87 and 97 km where substantial GW dissipation occurs. The cooling rate can exceed 30 K day −1 . In the lower region of wave dissipation, the strong downward sensible heat flux generates large heating. The net effects can be seen from the annual profile in Figure 5b. In general, the heat transport by GW dissipation induces considerable cooling throughout the mesopause region. This net cooling can significantly influence temperature structure and cools down the mesopause.

Heat Flux due to Waves of Different Scales
In the previous section, we have analyzed the transport induced by dissipating GWs over the entire spectrum that lidar data can resolve, including GWs with periods from a few minutes to 6 h and vertical wavelength from 2 to 15 km. Here we examine the contributions to fluxes from GWs in three different period intervals with the same data set. Based on linear GW theory, compared to low frequency GWs, high-frequency waves tend to propagate more vertically and air parcels oscillate more vertically. They are expected to contribute more to the sensible heat flux when experience dissipation but this is yet to be confirmed observationally. The three intervals we consider correspond to wave periods from 5 min to 1 h, 1-3 h, and 3-6 h. Vertical scale remains the same for each period interval in order to keep one dimension consistent. The method used to extract these three sets of perturbations is described in Section 2.2.
Shown in Figure 6 are the annual vertical fluxes calculated from three sets of GW perturbations with different periods. 〈w′θ′〉 is calculated from the same algorithm as seasonal flux except that the integral interval of frequency spectrum (ω 1 , ω 2 ) now depend on the period ranges of extracted wave perturbations. Profiles in each plot represents the annual fluxes averaged over 118 nights of lidar observations. As seen from medium, and high period waves, we conclude that shorter period GWs (high frequencies) contribute more to the heat transport and dissipate at higher altitude in the mesopause region.
In addition to the heat transport by GWs, the energy transport can be related to the wave total energy and energy flux velocity. The total energy density per unit volume is defined as where u = (u, v, w) is the velocity vector, u as the zonal wind, v as the meridional wind, 2 s c RT   is the speed of sound, and γ = c p /c v , c p and c v are the specific heats at constant pressure and volume, respectively. The first term on the right hand side is the wave kinetic energy; the second term represents the wave potential energy; the third term is the elastic potential energy and is much smaller than the previous two terms Walterscheid & Schubert, 1990). The total energy per unit mass is / With the wave vertical energy flux F z = 〈w′p′〉, the vertical component of energy flux velocity, which is equal to the wave vertical group velocity c gz in the mesopause region (Vincent, 1984;, can be calculated from: Among the 118 nights of lidar observations, 82 nights were operated with off-zenith directions, so zonal and meridional winds (u and v) were derived from LOS winds in the East and South directions, respectively. The horizontal wind perturbations (u′ and v′) were then derived with the same approach used for w′ and T′. (u′, v′, w′) were then applied to calculate the kinetic energy and T′ was used to estimate the potential energy. Summarized in Table 2 are the mean energy per unit mass E m , energy flux per unit mass / z F , and vertical group velocity c gz between 85 and 100 km for three different scale range waves. It is obvious that energy flux is larger for high-frequency waves and smaller for low frequency waves. The vertical group velocity can also be derived from the derivative of intrinsic frequency with respect to vertical wavenumber as (Fritts & Alexander, 2003): where k and l are horizontal wavenumbers, f is inertial frequency, and H is the scale height. Our results in Table 2 shows a good agreement with Equation 8 that c gz increases with frequency. The average c gz for high, medium, and low frequency waves are 4.15 ms −1 , 1.15 ms −1 , and 0.70 ms −1 , respectively. The group velocities are also consistent with previous studies (e.g., Chen et al., 2013;Lu et al., 2015a;Vincent, 1984).

Discussion and Conclusion
In this paper, we presented the calculation of vertical fluxes of sensible heat, enthalpy, and kinetic energy of dissipating GWs based on a Na lidar observation at ALO. This is the first time that the sensible heat flux in the mesopause region is measured in the Southern Hemisphere. Compared with the results reported in GL07, despite the opposite seasons, same feature appears in both locations, that sensible heat flux is maximum during solstices and minimum around equinoxes, but with differences in the altitudes and magnitudes of the maximum sensible heat flux. We have also analyzed the seasonal fluxes of enthalpy (〈w′θ′〉, or potential temperature flux) and energy flux (〈w′p′〉) for the first time. When the wavefield consists primarily of upward propagating GWs, the energy flux is positive and the enthalpy flux is either zero or negative depending on whether waves are dissipating. The sensible heat flux is also generally downward, except when there is very little dissipation and it becomes upward and proportional to the relative energy flux ( / w p p    ). Our results show that the enthalpy flux is downward throughout the mesopause region in all seasons, which is consistent with the well-known theoretical studies (Walterscheid, 1981;Weinstock, 1983). The downward heat transport generated by wave dissipation results in a significant cooling to the mesopause with a cooling rate as much as of 15 K day −1 .
We have also investigated the contribution of different scale GWs in downward heat transport. Results suggest that smaller period waves generate larger heat fluxes and dissipate at higher altitudes. The energy densities and energy flux are largest for long period waves. The upward energy flux decreases from 10 −2 W m −2 at 85 km to 10 −4 W m −2 at 100 km. The vertical group velocities are inferred from energy flux and total wave energy. The values are consistent with the knowledge that high-frequency waves propagate more vertically. The magnitudes of our inferred vertical group velocities are also comparable to those obtained with other techniques in the mesopause.
In most general circulation models, the GW parameterization primarily focuses on the wave drag through momentum flux. The thermal effect of dissipating GWs are not well represented. The measured heat flux and associated cooling rate, together with the vertical energy flux in this study can help us better understand the impacts of wave dissipation in the MLT and provide references for future wave parameterization in global models.
For energy flux measurement, one important assumption we have to make is that the measured frequency spectrum is a reasonable representation of the intrinsic frequency spectrum. There is generally some deviation from this because wav es are likely propagating opposite to the background wind, therefore having a intrinsic frequencies higher than their measured ground-based frequencies. This causes a shift of power toward lower frequency and an increase of spectral slope in observed frequency spectra compared with intrinsic frequency spectra. While the impact of this on the energy flux is estimated in this study, a more refined estimate can be obtained by using horizontal wind measurement and/or airglow imaging to determine the statistical distribution of wave propagation directions.

Appendix A: Sampling Error in Heat flux Calculation
〈w′T′〉 is a small quantity so careful analysis of its measurement errors is especially important. There are several potential sources of error in 〈w′T′〉 calculation: 1. Random errors ΔT and Δw due to photon noise contribute to the error in 〈w′T′〉 2. Statistical fluctuations of w′ and T′ cause fluctuations of 〈w′T′〉 and uncertainty of the mean value of 〈w′T′〉 3. Correlations between ΔT and Δw can create a bias in 〈w′T′〉 4. Systematic errors in T′ and w′ can also induce a bias in 〈w′T′〉 5. Limitation in temporal and spatial resolutions of w′ and T′ in resolving full GW oscillations can contribute to errors in 〈w′T′〉 Among these error sources, unbiased noises (1 and 2) can be estimated but not eliminated. They can only be reduced through averaging of more independent samples. Biases (3 and 4) can be corrected if they can be quantified. Error source 4 may cause both uncertainty as well as bias. Since w′ and T′ are nearly in quadrature and 〈w′T′〉 is close to zero for freely propagating GWs, accurately resolve small phase shifts from the quadrature may be needed for accurate 〈w′T′〉 measurements. Among these errors, one and two were addressed in Gardner et al. (2005) and Kudeki and Franke (1998). Three was studied by Su et al. (2008). The conclusion is that the error from two is dominant at typical lidar signal level, and requires tens of hours of averaging to reduce it down to an acceptable level, at which point the error from one and the bias from three becomes negligible. The bias in four comes from errors in retrieving temperatures and wind velocities due to insufficient knowledge of lidar system parameters. One such example is the correction of the nonlinear effect of photomultiplier tube. This was addressed recently by A. Z. , who showed a procedure to remove these effects objectively based on measurement. The data set used in this study have been corrected for this bias. Here we further examine the impact of error source 5.
At ALO, the typical lidar measurement mode is 60-s integration in zenith (Z), south (S) and east (E) directions sequentially as ZES, or as ZEZS. On some nights the lidar was pointed to zenith only. Vertical profiles of temperature and line-of-sight (LOS) wind velocities are derived at 500-m range resolutions. Only the zenith profiles are used for heat flux calculation. Thus the temperature and vertical wind measurements are available at a 3-min cadence in ZES mode, 2-min in ZEZS mode, and 1-min in zenith only mode. Since the period of GWs can be as short as the buoyancy period of about 5 min (ignoring the Doppler effect), this 3-min cadence in ZES mode is expect to have significant error for heat flux measurement of a high-frequency GW. On the other hand, we note that the goal here is not to measure the heat flux of an individual GW, but to obtain a statistically reliable long-term mean value. Over tens of hours, there are many GWs with varying periods, amplitudes and phases that contribute to the net 〈w′T′〉. It can be expected that the large error of heat flux measured for a single GW could be reduced when averaged over a large ensemble of random GWs.
This sampling issue is illustrated in Figure A1. For a 5-min wave, there are only two samples in one wave period at 3-min cadence. As shown in Figure A1a, the two samples in the first period happen to be at locations where positive w′T′ (both w′ and T′ are positive in the first sample and both are negative in the 2nd sample), resulting in a positive 〈w′T′〉. In the second period, both samples give negative w′T′ therefore a negative 〈w′T′〉 for this period. Heat flux calculated within each one wave period is very different from its true value, but the average over these two periods is much closer to the true value. For longer period waves, the sinusoidal variation is better sampled. As shown in Figure  We assume that the GW perturbations of temperature and vertical velocity are sinusoidal functions of time, as T′(t) = A T sin(ωt) and w′(t) = A w cos(ωt + ϕ). A small phase difference ϕ gives rise to a non-zero heat flux. We use τ to represent the integration time, which is the time the lidar spends to obtain one set of temperature and vertical velocity profiles; and Δt to represent the cadence, the time between the beginning time of two consecutive measurements. We also denote t 1 = t 0 + Δt as the time the first measurement starts, and t i = t 0 + iΔt as the starting time of the ith measurement. Each measurement is an average of the true quantity over τ duration. Thus the ith temperature and vertical velocity perturbations are 2 sin( ) sin sin , 2 2 The heat flux (HF) is the average of the above from N measurements is the true heat flux, and b and e represent two errors in HF that cause it to deviate from HF 0 , expressed as The error b is related to the finite integration time τ. Since 0 < b < 1, it causes a reduction in the HF magnitude. This reduction is smaller for smaller τ and smaller ω (longer wave period). This relationship is shown in Figure A2 for τ = 0.5 min, 1 min and 1.5 min. It shows that for the shortest 5-min wave, the magnitude reduction can exceed 10% with τ = 1 min but is < 1% for wave period longer than 20 min.
Error e is a sinusoidal function of t 0 with the amplitude max 1 sin( Δ ) .
When the sampling times are uniformly distributed (a reasonable assumption when the waves are randomly sampled many times), the mean value of e is zero so it is an unbiased error. e is not Gaussian-distributed and is more likely to be at the maximum value (either positive or negative) than near zero. Because of this, it is more appropriate to use e max than the rms to assess this sampling error. Generally, e max decreases as the number of samples N increases, because when the waves are sampled at more phase values the average of e is closer to the true value zero. In a special case when ωΔt = π, namely the wave period is twice the cadence, so e max does not decrease with the increase of N. However, in the real atmosphere, this condition can never be satisfied for an extended period of time so increasing N is always effective in decreasing e max . Figure A3 shows numerical calculated e max for two cadences of 1 min and 3 min, sampled over a single wave period for the period range from 5 min to 60 min. Thus the total samples N is larger for longer period waves because of longer sample duration. In the figure, e max is plotted as light gray lines. It reaches local maximums and minimums alternatively due to the nature of sinusoidal functions in Equation A8. These large fluctuations of e max however, are not realistic in the real atmosphere, because of natural fluctuations of wave periods. Therefore, we also plotted the running averages of e with a 2Δt window as red lines, which are more representative of the expected maximum error. When Δt = 3 min, this averaged e max varies from 35% for 5-min wave to 0.49% for 60-min wave. When Δt = 1 min, this averaged e max varies from 5.1% for 5-min wave to 0.31% for 60-min wave. This  shows that the cadence can contribute to a large error in HF for waves with period close to the cadence, but the error decreases quickly as the wave period increases. The larger e max for shorter-period waves is due to less sampling within one wave period.
More relevant to heat flux measurement is the errors for sampling over a fixed duration instead a wave period. As shown in Table 1, there are tens of hours of observation in each calendar month. Considering that GWs may not appear all the time, we chose 8 h as the sampling duration to estimate e max . This should give a conservative estimate of maximum possible error due to sampling. For a given sampling duration, shorter-period waves are sampled more than longer-period waves. This compensates the larger errors for shorter-period waves shown in Figure A3. The e max with 8-h sampling duration and Δt = 3 min is shown in Figure A4. The overall error is much smaller (except for the special case when ωΔt = π as discussed earlier), with the running mean |e| max at 0.39% for 5-min wave and 0.61% for 60-min wave, with the minimum of 0.16% at the period of 12.8 min.
The combination of the two errors b and e max gives the range of measured HF as max 0 HF HF 1 . sin This is illustrated in Figure A5 using b and e max shown with the red lines respectively in Figures A2 and A4 (8-h sampling, Δt = 3 min, τ = 1 min). Instead of plotting the percentage errors, we use nominal values A T = 5 K, A w = 1 ms −1 , and ϕ = 15° and show the heat flux with and without errors. In the figure, the black dashed line is HF 0 b, which shows a maximum reduction of the HF magnitude from 0.647 ms −1 to 0.566 ms −1 (12.5%) for the 5-min wave and near zero reduction (0.1%) for the 60-min wave. The red dashed lines are HF 0 (1 ± e max / sin ϕ). It shows that this error causes HF to vary between − 0.647 ± 0.0098 ms −1 (±1.5%) for the 5-min wave and − 0.647 ± 0.0154 ms −1 (±2.4%) with average of − 0.647 ± 0.0078 ms −1 (±1.2%) for the 60-min wave. Their combined effect as in Equation A10 are shown as the blue lines. Overall, the HF varies between 86.2% and 88.9% of the true value for the 5-min wave and between 97.5% and 102.3% for the 60-min wave, with an average of 97.7%-100.1%.
From the above error analysis, we can see that b is independent of the HF magnitude. e is dependent only on the phase different ϕ and independent of the amplitudes of T′ and w′. Overall the total fractional error is largely independent of the heat flux magnitude. As shown in Figure A5, the sampling error causes a very small decrease of HF magnitude, likely within 2% considering the lengths of data used in our calculation for all calendar months.
GUO AND LIU 10.1029/2020JD033825 14 of 16 Figure A5. Errors due to cadence Δt (black dashed line) and integration time τ (red dashed lines) and their combined error (blue lines). The black horizontal line is the true heat flux HF 0 . The error b due to finite integration time is shown as HF 0 b and the cadence error is shown as HF 0 (1 ± e max / sin ϕ). The value of b is the same as the red line in Figure A2 for τ = 1 min and the values of e max is the same as the red line in Figure A4 for Δt = 3 min and 8-h sampling duration.