Advanced hodograph-based analysis technique to derive gravity waves parameters from Lidar observations

An advanced hodograph-based analysis technique to derive gravity wave (GW) parameters from observations of temperature and winds is developed and presented as a step-by-step recipe with justification for every step in such an analysis. As a most adequate background removal technique the 2D-FFT is suggested. For an unbiased analysis of fluctuation whose amplitude grows with height exponentially we propose to apply a scaling function of the form exp(z/(ςH)), where H is scale height, z is altitude, and the constant ς can be derived by a linear fit to fluctuation profile and should be in a range 1–10. The 5 most essential part of the proposed analysis technique consist of fitting of cosines-waves to simultaneously measured profiles of zonal and meridional winds and temperature and subsequent hodograph analysis of these fitted waves. The linear wave theory applied in this analysis is extended by introducing a wave packet envelop term exp(−(z− z0)/2σ) that accounts for limited extent of GWs in observational data set. The novelty of our approach is that its robustness ultimately allows for automation of the hodograph analysis and resolves many more GWs than it can be inferred by manually applied hodograph technique. 10 This technique allows to unambiguously identify upand downward propagating GW and their parameters. This technique is applied to unique lidar measurements of temperature and horizontal winds measured in an altitude range of 30 to 70 km.


Introduction
It is generally accepted that atmospheric gravity waves (GW) produce global effects on the atmospheric circulation from the surface up to the mesosphere and lower thermosphere (MLT) region (e.g., Fritts and Alexander, 2003;Alexander et al., 2010;Becker, 2017).Well known tropospheric sources for these waves are the orography (flows over mountains), convection, and jet imbalance (e.g., Subba Reddy et al., 2005;Alexander et al., 2010;Mehta et al., 2017).When propagating upwards, GW dissipate and thereby deposit their momentum starting from the troposphere and all the way up to the MLT.This process is referred to as GW-forcing and plays a key role in the global circulation.The most of the climate models are not able to resolve these small-scale waves (i.e., waves with horizontal wavelengths typically shorter than 1000 km) (e.g., Kim et al., 2003;Geller et al., 2013).That is why these waves and their dissipation (and also their interaction with each-other and with the background flow) are often called "sub-grid scale processes" (e.g., Shaw and Shepherd, 2009;Lott and Millet, 2009).In order to account for the influence of GW these models need to rely on various parametrizations.To construct a proper parametrization one has to describe GW frequencies, wavelengths, and momentum flux over the model coverage zone (e.g., Alexander et al., 2010;Bölöni et al., 2016).
Our knowledge about gravity wave parameters can be improved by means of high resolution measurements of atmospheric GW.Ideally, the measurement range should cover the entire path of the waves starting from their sources in the troposphere to the level of their dissipation, that is up to the MLT region.Such type of measurements ultimately faces high experimental challenges which explains why we still do not have satisfactory and conclusive observational data on these processes.
In the altitude range of the mesosphere only few observation techniques exist.In the last decades the only source of highresolution GW observations based on both temperatures and winds in the stratosphere and mesosphere region were rocket soundings (see e.g., Schmidlin, 1984;Eckermann and Vincent, 1989;Lübken, 1999;Rapp et al., 2002, and references therein).
Rocket measurements with e.g., falling spheres can provide vertical profiles of horizontal winds and atmospheric temperatures and densities with altitude resolution of about 1-10 km.
Satellite-borne remote sensing techniques can provide excellent global coverage, their observations deliver unique horizontal information about GWs (see e.g., Alexander et al., 2010;Alexander, 2015;Ern et al., 2018), but they base solely on temperature observations.Ground-based radar systems are able to measure winds at heights 0-30 km and 60-100 km.From the altitudes between 30 and 60 km radars do not receive sufficient backscatter and, therefore cannot provide wind measurements in this region.
While the vertical wave structure can be resolved from rocket profiles, the long and irregular time intervals between successive launches prevent the study of temporal gravity-wave fluctuations over a larger time span (Eckermann et al., 1995;Goldberg et al., 2004).
Recent developments in lidar technology give us new possibilities to study GW experimentally on a more or less regular basis and resolve spatial sales of 150 m in vertical and temporal scales of 5 min (e.g.Chanin and Hauchecorne, 1981).In particular the day-light lidar capabilities allow for long duration wave observations (e.g., Baumgarten et al., 2015;Baumgarten et al., 2018).The new Doppler Rayleigh Iodine Spectrometer (DoRIS) additionally to the established lidar temperature measurements yields simultaneous, common volume measurements of winds (Baumgarten, 2010;Lübken et al., 2016).This combination of capabilities makes lidar data unique.
All those quantities, i.e. winds and temperature, when measured with high temporal and spatial resolution, reveal structuring at scales down to minutes and hundreds of meters.In our analysis technique we aim solely at such fluctuations which are generated by GW.By applying a proper data analysis technique one can extract several important parameters of GW from the advanced lidar measurements.
In this paper we describe a newly developed analysis technique which allows for derivation of GW parameters such as vertical wavelength, direction of propagation, phase speed, kinetic and potential energy and momentum flux from the advanced lidar measurements.We aim at presenting a step-by-step recipe with justification of every step in such an analysis.Every single steps if considered independently, are in general well known.The strength and novelty of our work is their combination and some justification on their importance and how they affect analysis results.The paper is structured as follows.In the next section a short description of lidar measurement technique is given.Section 4 describes the new methodology in detail.Finally, in section 5 geophysically meaningful quantities are deduced from the analyzed data which also demonstrates the capabilities of the introduced analysis technique.Theoretical basis used by the data analysis technique is shortly summarized in section 3 and extended in Appendix A.

Instrumentation
The ALOMAR Rayleigh-Mie-Raman lidar in northern Norway (69.3 • N, 16.0 • E) is a Doppler lidar that allows for simultaneous temperature and wind measurements in the altitude range of about 30 to 80 km.The lidar is based on two separate pulsed lasers and two telescopes (von Zahn et al., 2000).Measurements are performed simultaneously in two different directions, typically 20 degrees off-zenith towards the North and the East by pointing the telescopes and the outgoing laser pulses in this direction.
The diameter of each telescope is about 1.8 m and the average power of each laser is ∼14 W at the wavelength of 532 nm.
Both pulsed lasers operate with a repetition rate of 30 Hz and are injection seeded by one single CW-laser that is locked to an Iodine absorption line.The light received by both telescopes is coupled alternatingly into one single polychromatic detection system.Temperatures and winds are derived using the Doppler Rayleigh Iodine Spectrometer (Baumgarten, 2010).As the measurements discussed below are performed also under daytime conditions we process the data as described in Baumgarten et al. (2015).Measurements by the lidar were extensively compared to other instruments showing the good performance of the lidar system (Hildebrand et al., 2012;Lübken et al., 2016;Hildebrand et al., 2017;Rüfenacht et al., 2018).The lidar data are recorded with an integration time of 30 seconds and a range resolution of 50 m.The data are then integrated to a resolution of 5 minutes and 150 m and afterwards smoothed with a Gaussian window with a full width at half maximum of 15 minutes and 0.5 km is performed.For calculation of horizontal winds from the measured line-of-sight winds we assume that the vertical wind component is equal to zero.Importantly, the estimated uncertainty imposed by this assumption is negligible and does not affect final results of our analysis.The hydrostatic temperature calculations were seeded using measurements from the IAP mobile Fe resonance lidar and the temperatures from both lidar systems were then combined by calculating an error weighted mean (Lautenbach and Höffner, 2004).

Brief theoretical basis
A GW-field consists of various waves with different characteristics.An attempt to describe this system as a whole is made, for example, by Stokes analysis (e.g., Vincent and Fritts, 1987;Eckermann, 1996).In this work we do not try to describe bulk fluctuations, but rather to extract the single most dominant quasi monochromatic (QM) gravity waves (GW) from the set of the observed fluctuations.The advantage of this approach is that it allows us to describe these selected waves as precisely as possible by the linear theory of GW.Moreover, the main idea of our retrieval is to find GW-packets where fluctuations of the both wind components, i.e. zonal and meridional wind (u and v ) as well as temperature fluctuations T show the same characteristics, i.e., belong to the same wave-packet.This requirement ensures, that our analysis only accounts for wave structures and not for those created by accompanying dynamical processes like turbulence or other wave-like structures created by e.g., temperature inversion layers (e.g., Szewczyk et al., 2013).
For this analysis we use the assumption, that a wave packet at a fixed time point and in a limited altitude range can be considered as QMGW, i.e. dispersion within one wave packet is neglected.Also we assume, that all the observed parameters (T, u, v) reveal fluctuations (T , u , v ) at the same wavelength.Mathematically it can be written in the following form (see Appendix A for more details): where ϑ refers to either of temperature (ϑ ≡ T ), zonal or meridional wind components (ϑ ≡ u or ϑ ≡ v); prime variables describe fluctuations (T , u , v ) and ϑ is amplitude of those fluctuations; ϕ ϑ is phase shift; m is the vertical wave number (λ z = 2π/m is vertical wavelength) and H is the scale height.
Eq. 1 is an ansatz which describes an ideal monochromatic GW under the conditions of conservative propagation in a constant background.Similar description of GW propagation is widely used in the literature (see e.g., Gavrilov et al., 1996).
Since most GW propagate oblique through the field of view of the ground-based instruments, they appear in the observations as waves of a limited vertical extent, i.e. as wave packets.Although, any and also vertically propagating waves might appear in the nature in the form of wave packets rather than continuous wave of quasi-infinite length.Therefore we extend the Eq. 1 by introducing a wave packet envelop term exp(−(z − z 0 ) 2 /2σ 2 ) that accounts for limited presence of the GW-packet in observations: where σ is a factor that describes width of wave packet and z 0 is the altitude of maximum of wave envelope (its central altitude).
Following e.g., Cot and Barat (1986) or Gavrilov et al. (1996), the horizontal propagation angle of QMGW can be defined as: where ξ is the azimuth angle of wave propagation direction and When v > u, n = 0 and 2 for F uv > 0 and F uv < 0, respectively.This implies, that for ϕ u − ϕ v = π/2 propagation direction can be 0 or 180 degrees, i.e. northward or southward if v > u and eastward or westward if v < u.The sign of m in Eq. 1 and 2 shows the vertical propagation direction: m < 0 for upward and m > 0 for downward propagating GW.
This theoretical basis allows to describe the main GW-parameters and to derive them from observations.However in practice, noisy data and/or insufficient resolution of measurements may lead to large uncertainties when applying these equations directly to the measured time series.Therefore, most common technique, based on linear theory of gravity waves to derive propagation direction, intrinsic frequency and phase velocity of GW from ground-based observations is the hodograph method (e.g., Sawyer, 1961;Cot and Barat, 1986;Wang and Geller, 2003;Zhang et al., 2004;Baumgarten et al., 2015).The hodograph technique explicitly utilizes the following polarization relations of GW for winds and temperature.
For mid-and low-frequency GW the velocity perturbations in propagation direction and perpendicular to this direction are related by the polarization relation (e.g, Gavrilov et al., 1996;Fritts and Alexander, 2003;Holton, 2004): where v ⊥ is complex amplitude of wind fluctuations in the direction perpendicular to the direction of propagation and u is amplitude of wind fluctuations along the propagation direction.f = 2ΩsinΦ is Coriolis parameter and ω is intrinsic frequency.
That is, for a zonally propagating wave v ⊥ is the meridional velocity amplitude.
If we assume, that θ/θ = T /T 0 (Fritts and Rastogi, 1985;Eckermann et al., 1998), the temperature amplitude is related to the parallel wind amplitude for a wave propagating in zonal direction as (e.g, Hu et al., 2002;Geller and Gong, 2010): where k h = 2π/λ h is the horizontal wave number and λ h is the horizontal wavelength of the QMGW; θ/θ are potential temperature perturbations; T 0 and g are the background temperature and the acceleration due to gravity averaged over the altitude range of the QMGW.
To summarize, the basic theory briefly described in this section allows to derive the main GW parameters: intrinsic frequency, amplitude and direction of propagation.From these parameters one can derive a more extended set of GW parameters as summarized in App. A.

Retrieval algorithm
In this section we describe the procedure to derive wave parameters from the measured lidar data.For our analysis we need simultaneously measured wind and temperature profiles.Technically we can extract wave parameters from a single measurement, that is using two wind and one temperature profile.However, for a robust estimation of the atmospheric background we need a several hours long observational data set.

Separation of GW and background
The first step is to remove the background from the measured data.The background removal procedure plays a key role in GW-analysis techniques and may even lead to strongly biased results.The main reason for this is that the most analysis techniques rely on fluctuation's amplitudes remaining after subtraction of background to infer wave energy (e.g., Rauthe et al., 2008;Ehard et al., 2015;Baumgarten et al., 2017;Cai et al., 2017, and others).Since GW energy is proportional to amplitude squared, any uncertainty in the background definition ultimately leads to large biases in estimation of GW-energy.
We define the background as fluctuations with periods and vertical wavelengths longer that typical GW parameters.This means, that tidal fluctuations and planetary waves are attributed to background.Tides periods that are integer fractions of a solar day.Semidiurnal tides have period of 12 hours and coriolis period (2π/f ) at 69N is 12, 8 hours.Thus, only doppler shifted GW can reveal periods longer than ∼ 12 hours.From other site, typical vertical wavelengths of GWs was summarized in Table 2 of Chane-Ming et al. (2000) and not exceed 17 km.Thus, we define the background as wind or temperature fluctuations with periods longer than 12 hours and vertical wavelengths longer than 15 km.To extract such a background from measurements we apply a low-pass filter to the altitude vs time data.Specifically, we use the two dimensional fast Fourier transform (2D-FFT) (e.g., González and Woods, 2002) and, after blocking the specified high frequencies and short wavelengths, and applying the inverse 2D-FFT, we finally construct the background.Advantage of this method is that it simultaneously accounts for both variability in space and time.After subtracting the derived background from the original measurements we obtain the wind and temperature fluctuations which have periods shorter than 12 hours or wavelengths sorter than 15 km and supposedly produced by gravity waves.This procedure is demonstrated in Fig. 2, 3, and 4 for temperature, zonal, and meridional wind, respectively.The upper, middle, and lower panels represent the original measured quantities, estimated background, and the resultant fluctuations, respectively.These time-altitude plots consist of many single-time ("instant") altitude-profiles which are further analyzed individually.More specifically, fluctuations T , u , and v are analyzed with our automated hodograph method.
We also performed a robustness test to check how different background removals influence our advanced hodograph-based method.To derive the background (for both wind and temperature data) we made use of (a) running mean with different smoothing window lengths, (b) different splines, and (c) constant values in time.It turned out that our analysis results were near identical for all these different backgrounds.The new technique is not sensitive to the background derivation schemes and may even allow to skip this step from the analysis.A more in depth analysis showed, that the robustness to the background removal is a consequence of the analysis approach.We only search for waves which are prominent simultaneously in temperature and both wind components.Even though we are confident in the robustness of our technique to the various background derivation methods, we consider the 2D-FFT based approach as the one most adequate for this purpose.

Scaling of fluctuations
Under the assumptions of conservative propagation (i.e., without wave breaking and dissipations) and a constant background the amplitude of fluctuations increases with altitude as exp(z/(2H)).In the real observations, since waves cannot freely propagate throughout the atmosphere, the amplitude of the fluctuations increases with altitude as exp(z/(ςH)), where the coefficient ς ≥ 2 is derived from the observed data.The exponential growth, however also affects any analysis, in particular wavelet analysis, since normalization is always applied.The growing amplitude works as a weighting function and, therefore, the largest amplitudes will dominate the analysis (e.g., spectrum), thereby hiding the small-amplitude waves (see also Wright et al., 2017, who pointed out to similar effect in the satellite data).This effect, in particular, prohibits analysis of small-scale features at lower altitudes.Scaling the fluctuations by 1.0/ exp(z/(ςH)) yields fluctuations with comparable amplitudes over the whole altitude range.For the observations presented here we derived ς = 2.15.Note however, if further analysis requires treatment of fluctuation amplitudes, this scaling must be either taken somehow into account (e.g., by appropriate normalization) or removed (by applying inverse scaling) as we do in Sec.4.7.

Detection of wave packets
Starting from this point we only analyze the altitude-profiles at every time step.At every time step we have measured profiles of wind and temperature which are split in fluctuations and background profiles.
First, we search for dominant waves in both altitude and wave number domains.For this purpose we apply the continuous wavelet transform (CWT) to every profile of the extracted fluctuations.We use a Morlet wavelet of the sixth order (Torrence and Compo, 1998) and apply it to vertical profiles of wind and temperature fluctuations.Similar procedure was also applied by Zink and Vincent (2001) and Murphy et al. (2014).By applying wavelet analysis they define regions from which the Stokes analysis (e.g., Vincent and Fritts, 1987;Eckermann, 1996) is further evaluated with a better precision.We note here, that their results rely on accuracy of wavelet transform and on assumption that wave signatures are well separated from each other and clearly resolved by CWT.
An example for the resulting scalograms of one time step is shown in Fig. 5.These scalograms are normalized to unity to make spectral signatures comparable between the different fluctuations.In zonal wind and temperature fluctuations a clear peak between ∼40 and ∼55 km with a vertical wavelength of approximately 10 to 15 km can be seen.Both wind components reveal peaks below ∼40 and above ∼60 km with wavelengths of about 5 km.As a next step we combine these wavelet spectra and construct a single scalogram that reflects the features common for all three components.We calculate the product of all three spectra and define this as the combined spectrum.Note, that Zink and Vincent (2001) and Murphy et al. (2014) used sum of scalograms of both wind components.The combined scalogram in Fig. 6 reveals one large (around 10 km wavelength) and two smaller (near 35 and 70 km altitude) regions with weaker wave amplitudes.The larger region is relatively broad and reveals a vertical wavelength increase with increasing altitude.This can be due to two reasons: 1) it is one wave packet with changing vertical wavelength due to variable background or 2) it is a sum of two wave packets with overlap at around 50 km altitude.This uncertainty is difficult to resolve just using information from wavelet transform.To resolve this ambiguity we developed a sequence of further analysis steps and only use these CWT results as an input (zero guess) for further analysis.

Fitting of linear wave theory
In this step we fit a wave-function to all three measured profiles, i.e. u (z), v (z), and T (z).The wave-function was derived from the linear wave theory as summarized in Sec. 3 and App. A. Note, that wave-function described by Eq. 1 includes scaling factor exp(z/(2H)).After applying the step of our algorithm described in Sec.4.2 (scaling), we get rid of exponential growth in fluctuations profiles and, thereby exclude this factor from wave equation.Therefore, as wave-function that we fit to the profiles u (z), v (z), and T (z) we only use the remaining part: where ϑ refers to u, v and T .
The fit can be performed using a least square regression algorithm implemented in numerous routines.The data (measurements) to which the wave-function is to be fitted are the three profiles u (z), v (z), and T (z).The fit must converge for all three profiles to be qualified as successful.The free fitting parameters are central altitude of the wave packet z 0 , width of the wave packet σ, wavelength of oscillations in this wave packet λ z , amplitude of fluctuations in the wave packet | ϑ|, and the phase shift ϕ ϑ .Initial guess for parameters λ z , z 0 , and σ is estimated from the wavelet scalogram derived in the previous step.
Zero guess for the amplitude | ϑ| is directly derived from the fluctuation profile as maximum amplitude in the height range z 0 ± λ z /2.Initial value for the phase shift ϕ ϑ is taken randomly.
Thus, to derive first set of initial paramteres λ z , z 0 , and σ we start with the larger area encircled by the dashed lines in Fig. 6 and and pick up the values λ z =12 km, z 0 =45 km, and σ=15 km.The initial amplitude for e.g., zonal wind fluctuations estimated from the red profile in Fig. 7a | u| =10 ms −1 .The fit of Eq. 1 to the temperature and two wind profiles will yield set of parameters that describe a wave packet: Thus, the updated values for this demonstration case are z 0 = 49 km, λ z = 11 km.
We recall that the introduced in Sec. 3 vertical extend of wave packet exp(−(z − z 0 ) 2 /2σ 2 ) is essential for analysis of observations which cover an altitude range of approx.50 km and thus much longer than a wavelength and the expected scale of amplitude variations.
Similar way of deriving initial guess parameters was, for example, implemented by Hu et al. (2002), who used power spectrum to define dominant waves.However, since their observations only cover 20 km altitude, they do not need to consider thickens of wave packet.Hu et al. (2002) simply assumed, that wave packet covers the entire altitude range of their observations.
Obviously, such an assumption is not valid if observations cover an extend altitude range like in our study.
Step 4.3 and, in particular Fig. 6, clearly support this statement.
Generally speaking, intrinsic frequency and propagation direction can be estimated from the obtained fitting results by applying Eq. 3 and A5 respectively.However, by testing different simulated and measured data we concluded that for GW with intrinsic periods larger than ∼1 hour the hodograph analysis yields more accurate results than those based on the fitting of Eq. 6.Therefore, the described fitting procedure is only used to precisely derive the altitude z 0 and the vertical wavelength λ z = 2π/m of the wave packet, which are smeared in the spectrogram (Fig. 6).Thus, we continue analysis using the hodograph technique.

Hodograph method
According to the theory (see App.A and e.g., Sawyer, 1961;Cot and Barat, 1986;Wang and Geller, 2003;Zhang et al., 2004;Baumgarten et al., 2015) the u and v fluctuations form an ellipse if the intrinsic period is between ∼1 and ∼12 hours.For Higher frequency GW, i.e. with periods below ∼1 hour the fluctuations form a line, as the influence of the Coriolis force is negligible.For low frequency GW, i.e. those with periods close to the Coriolis period (2π/f ) the fluctuations reveal a circle.
To extract the essential parameters of the wave packet found in previous steps we apply the hodograph analysis around the center of the wave packet (e.g.Baumgarten et al., 2015).Fig. 1 schematically illustrates this method.In the center of the wave packet the QMGW produces fluctuations in zonal and meridional wind components with equal vertical wavelengths but different phases and amplitudes, which is described by Eq. 1.The left panel of Fig. 1 shows the u and v wind fluctuations as a function of altitude.One can see several oscillation periods centered around ∼40 km altitude.If we select one full wave period around the center altitude of the QMGW, i.e. from z 0 − λ z /2 to z 0 + λ z /2, and plot u versus v , we get an ellipse as shown in the right panel of Fig. 1.The selected height range with one wave period is marked in Fig. 1a by the shaded area.The major axis of the ellipse is oriented along the wave propagation direction.
In order to minimize an error in the hodograph analysis due to presence of other waves (Zhang et al., 2004), we apply a vertical band-pass filter to all three profiles and thereby remove waves with wavelengths shorter than λ z /2 and longer than 2λ z .For example, if the vertical wavelength obtained in the previous step is 10 km, we remove waves with wavelengths shorter (longer) than 5 (20) km.Using a two dimensional least square fitting software we find the best fit parameters that satisfy the ellipse equation (Fitzgibbon et al., 1996).The fitting procedure is sensitive to the data quality and if for example the data is far away from an elliptical shape, the fitting procedure does not converge.Only if the ellipse was successfully fitted, we extract further wave characteristics from this data set.
The vertical propagation direction of the wave is unambiguously determined by the rotation direction of the zonal wind versus meridional wind hodograph.In the northern hemisphere the (anti-) clockwise rotation of the hodograph indicates a (downward) upward propagating wave.
The rotation direction of the hodograph is defined as a phase angle change of either u or v from the bottom level to the top level over the height region An additional hodograph of the parallel wind fluctuations versus temperature fluctuations is used to resolve an ambiguity in horizontal propagation direction that arises from the orientation of the ellipse in Fig. 1b.

Optimization of results
If in the previous step the rotation of hodograph does not make a full 360 • cycle, this suggests either an inconsistency in the hodograph results (sec.4.5) and the wave fit (sec.4.4) or the vertical extent of the wave packet is smaller than its vertical wavelength.In such a case we apply a correction to the vertical wavelength derived in the step 4.4.This correction to the vertical wavelength is found by forcing the hodograph to close the full 360 • cycle and calculating the additional vertical length resulted from this extra rotation.If the new (corrected) wavelength λ z differs significantly from λ z obtained before (sec.4.4), we repeat the hodograph analysis using the new (corrected) wavelength.

Calculation of GW parameters
The ratio of the major and minor ellipse axes is further used to derive the intrinsic frequency of GW (App.A).The analysis presented so far allows to derive the intrinsic wave frequency, vertical wavelength, and up-or downward propagation.The horizontal direction of propagation is along the major axis of the ellipse with a remaining uncertainty of 180 • .To resolve this ambiguity we further use temperature fluctuations profile as described in sec.3. Specifically, we construct a hodograph from temperature fluctuations and wind fluctuations along the wave propagation direction, i.e. parallel to the derived wave vector.
The rotation direction of this new hodograph finally defines the direction of the wave vector: for upward propagating GW clockwise or counterclockwise rotation indicates eastward or westward direction, respectively.For downward propagating GW opposite direction has to be used (Hu et al., 2002).
Knowing all these wave parameters and applying linear wave theory we derive further wave characteristics as described in App. A. Note, that as mentioned in Sec.4.2, at this point the fluctuation amplitudes must be rescaled back to their original growth rate with altitude using the derived scaling parameter ς, to legitimate their use for e.g., estimation of wave energy.
Fig. 7(a-c) show the fluctuations and two hodographs defined from the two maxima shown in Fig. 6. Results obtained from this example are summarizes in Table 1 (first column).

Iteration process
After the first QMGW is identified in all three profiles, it is subtracted from those data.We repeat the procedure described above for all of the maxima seen in the combined spectrogram (Fig. 6).In order to avoid over fitting, we limit our analysis to maximum of 20 waves per one time step.As it will be demonstrated in the next section that we never reach this limit.
In the given example 5 waves were detected.The first 3 waves are demonstrated in Fig. 7 and the obtained parameters are summarizes in Table 1.In this example we found that a wave with a vertical wavelength of 16.4 km propagates upward and against background wind in the altitude range from 56 to 73 km.In the altitude range of 44 to 54 km a wave with 11 km vertical wavelength propagates downward and with nearly the same direction as background wind.The analysis indicates that the broad maximum in the combined spectrogram (Fig. 6) was produced by the sum of two wave packets with different characteristics.

Reconstruction of 2D fields
Finally, this algorithm for a single point in time is subsequently applied to all time points of the entire data set shown in Fig. 2, In this section we demonstrate on a real data set how our analysis works and results are summarized in form of different statistics.
The data used in this study were obtained from 09 to 12 January 2016.During this time period a strong jet with wind speeds of more than 100 m/s was observed at an altitude range of 45 to 55 km (Fig. 3 and Fig. 4).During this period maps of the horizontal winds extracted from ECMWF-IFS (European Centre for Medium-Range Weather Forecasts -Integrated Forecasting System) showed a strong polar Vortex with wind speeds of more than 160 m/s at the vortex edge.The Vortex was elongated towards Canada and Siberia and its center displaced towards Europe.ALOMAR was located roughly below the Vortex Edge where the Polar Night Jet was located south of ALOMAR, at about 60°N with wind speeds of more than 160 m/s.
After applying the new analysis technique to the ∼60 hours measurements shown in Figs. 2, 3, and 4 we obtain the following results.
The number of detected waves per altitude profile is summarized in histogram Fig. 8.In 645 out of 715 altitude profiles we find at least one height range with a dominant GW where the hodograph analysis provides a reliable result.We recall that the analysis technique allows for up to 20 waves in a single profile.The total number of the detected waves amounts to 4507.It is seen that the majority of profiles yields 5 to 10 waves and none of them reaches the 20-waves limit.From the rotation direction of the velocity hodographs we derive that 32.3 % of all the detected waves propagate downwards.
Fig. 9 shows details of the wave packets as functions of altitude and separated for upward and downward propagating GW.
First plot shows the number of wave center altitudes (z 0 ) and does not consider the vertical extent of the wave packets (vertical wavelengths).The latter is taken into account in the middle panel which shows the mean fraction of the profile where a wave packet is present (any part of the wave, center or tail).We find that the most active regions (in terms of number of GW) are ∼32 to 40 km and 58-64 km.The altitude region between ∼40 and 55 km contains the smallest number of the detected waves.
It is interesting to compare these results with mean background wind shown in the rightmost panel of Fig. 9.It is obvious that the minimum in the wave activity as deduced by our analysis technique is co-located with the maximum of mean zonal wind as well as the background temperature.
The existence of downward propagating waves was reported earlier from observations by different methods (e.g., Hirota and Niki, 1985;Gavrilov et al., 1996;Wang et al., 2005) and also from model simulations (e.g., Holton and Alexander, 1999;Becker and Vadas, 2018).However, reported amount of downward propagating GW is very variable since most of observations were done at different altitudes or latitudes.Hu et al. (2002) found 223 (71 %) waves propagating upwards and 91 (29 %) downwards in the altitude range 84-104 km, which is in accord with our results.Gavrilov et al. (1996) reported that up to 50 % of the detected waves propagate downwards in the altitude range 70 to 80 km.In the troposphere and lower stratosphere (below 20 km) Sato (1994) reported less than 10 % downward propagating GW and Mihalikova et al. (2016) reported 18.4 % during wintertime and 10.7 % during summertime.From rocket observations of zonal and meridional wind components with a vertical resolution of 1 km in altitude range 30 to 60 km Hirota and Niki (1985) found in middle and high latitudes about 20 % of downward propagating GW, and 30-40 % in low latitudes at northern hemisphere stations.At the only southern hemisphere station (Ascension Island) 36 % of downward propagating GW were observed (Hirota and Niki, 1985).Hamilton (1991) found from rocketsonde observations of wind and temperature in the 28-57 km height range at 12 stations (spanning 8 • S to 76 • N) different fractions of downward propagating GW spanning from 2 % to 46 % depending on latitude and season.Wang et al. (2005) reported that approximately 50 % of the tropospheric gravity waves show upward energy propagation, whereas there is about 75 % upward energy propagation in the lower stratosphere.From their radiosonde observations authors demonstrate that the lower-stratospheric fraction of upward energy propagation is generally smaller in winter than in summer, especially at midand high latitudes.Thus, our finding of 32.3 % downward propagating GW reasonably agrees with other experimental data.
We note, that the observed downward or upward propagating GW are instantaneous observations, which means that we have no information about the fate of the observed waves.I.e., we cannot estimate the percentage of waves which ultimately get to the ground.
To investigate the time and altitude dependence of the GW detected by our hodograph technique, we reconstructed the temperature and the wind fluctuation fields from the derived waves parameters using Eq. 1. Fig. 10 shows the result of this reconstruction for the temperature fluctuations separated for upward and downward propagating GW.Contour lines show the background zonal wind velocity.We recall that the analysis technique treats every single altitude-profile independently and, therefore the influence of neighboring profiles is only due to time averaging.It is therefore remarkable, that the joint field of reconstructed GW shown in Fig. 10 builds up a consistent picture.Thus one can recognize for instance, waves packets of several hours duration.In some cases phase lines of waves follow the background wind.For example on 11 January after 18:00 UT at altitudes between 54 and 63 km a maximum of temperature fluctuations of upward propagating waves follows the contour line of a zonal wind of 60 m/s.
We use similar representations to investigate the temporal variability of any other of the derived GW properties.For example Fig. 11 summarizes the obtained intrinsic periods of GW throughout the measurement.On the one hand these figures demonstrate high variability, but on the other hand they also show regions of consistent picture.For example, on 11 January after 21 UT at altitudes between 54 and 62 km one can see wave period of about 7 hours for ∼ 2 hours.The analysis allows studying the temporal and altitude variation of the wave periods, e.g.upward propagating low period waves with large vertical wavelengths are often found above the jet maximum.
In Fig. 12 we show distributions of the derived GW parameters for all identified waves.One remarkable feature seen in these histograms is that the distributions of wavelengths and phase velocities reveal very similar shapes for up-and downward propagating waves.The distributions of intrinsic periods show quite different shapes, i.e. dissimilar kurtosis and skewness, for up-and downward propagating GW.These histograms also demonstrate limitations of the presented analysis.Only a few waves with intrinsic periods smaller than 1 h or with vertical wavelengths below 1 km are detected.This is likely caused by the smoothing of the lidar data with a Gaussian window of 15 minutes and 0.5 km rather than to the hodograph method itself.Waves with vertical wavelengths above ∼15 km were likely associated with background fluctuations when applying the 2D-FFT.
The distribution of phase velocities in Fig. 12 demonstrates that the velocities are below 60 m/s with a maximum of occurrence at ∼10 m/s.Matsuda et al. (2014) estimated horizontal GW phase velocities from airglow images.Their waves had periods below ∼1 hour and revealed phase speeds between 0 and 150 m/s.Among those waves, ∼70 % showed phase speeds between 0 and 60 m/s.In our case, the observed waves periods have maximum in the range 4 to 5 hours and, as expected from Eq. A10, the horizontal phase velocity is also lower than those, reported by Matsuda et al. (2014).
Another way to check the consistency of our technique is to look at the spectrum of fluctuations before and after analysis.As an example, Fig. 13 shows Fourier spectra of the temperature fluctuations calculated in time domain.The measurements and analysis results are represented by blue and orange lines, respectively.We recall that the analysis is made in spatial domain, that is it only deals with altitude profiles of fluctuations.Close similarity in both spectra which were calculated in time domain, that is across the analyzed profiles, suggests that the reconstructed two dimensional (time vs altitude) GW-field does not significantly deviate from the observed one.The reconstructed field indeed reflects the main GW-content and, therefore, in this respect it may be qualified as lossless algorithm.
Next, we analyze and sum up the wave energetics.Fig. 14 shows the derived kinetic (E kin ) and potential (E pot ) energy densities as well as their statistical basis.The altitude dependence of the energy distribution is shown as color coded two dimensional histogram (color bar on the right hand side defines number of waves).For these histograms all waves were used, i.e. both propagating up-and downwards.Fig. 14 also shows the mean energies separated for up and downward propagating waves.We find that E kin of downward propagating GW is lower than E kin of upward propagating waves and E pot is nearly identical for up-and downward propagating GW.The standard method to derive E kin and E pot from ground based observations is to average bulk wind and temperature fluctuations and apply Equations A13 and A15.Fig. 14 shows that the fluctuationbased method reveals a good agreement with mean profiles derived from our new retrievals.Obtained results for averaged E kin and E pot are also in agreement with mean winter profiles measured at ALOMAR observatory, summarized by Hildebrand et al. (2017).
The directions of background wind and wave propagation are summarized in Fig. 15 as polar histograms.The leftmost part, i.e.Fig. 15a shows a histogram of the background wind at the time/altitude of every hodograph.The analysis shows that in almost all cases the wind in the vicinity of detected waves blows towards the east-north-east with a mean speed of about 70 m/s.
Figs. 15b and 15c show polar histograms of the detected upward and downward propagating waves.From the color code we see, that the horizontal phase speed of the upward propagating waves is in general larger than that one of the downward propagating waves.Downward propagating waves reveal a rather uniform spatial distribution whereas upward propagating waves prefer to propagate against the background wind.
To address the question at which vertical angles the GW propagate we show histograms of the angle between the group velocity vector and the horizon (β; Eq.A12) separated for up-and downward propagating waves in Fig. 16.In the beginning of this section we noted that our analysis reveals that ∼30 % of all the detected waves propagate downwards.From the histograms we note that this difference of up-and downward propagating waves is mostly due to waves propagating at shallow angles of less than ∼1 degree.GW with larger vertical angles are found in same numbers for upward and downward propagating waves.
The vertical group velocities c gz estimated using Eq.A11 are summarized in Fig. 17 for up-and downward propagating waves.Since the vertical group velocity depends on wave periods, we split the histograms in two groups of longer and shorter than 8 hours.These results show for instance, that all low frequency waves (in the range of frequencies considered in our study) reveal small vertical group velocities.Waves with periods shorter than 8 hours show a somewhat more complicated picture.The vertical group velocities of downward propagating waves exceed those of upward propagating GW for waves that propagate in the direction of the background wind.In turn, upward propagating GW reveal highest vertical group velocities if they propagate against the background wind.The vertical group velocities c gz are at least two times lower than vertical phase speeds (c z , not shown here).The values of the vertical group velocity imply, that if waves propagate from the ground to the altitude where they were observed, they need 6 to 14 (2 to 4) days if they have period longer (shorter) than 8 hours.Somewhat similar time scales for GW to reach the lower stratosphere were reported by Sato et al. (1997) whose group velocity for waves with period of 17 h were 1.7 km/day.
Finally we show vertical fluxes of horizontal momentum (see Eq. A18) averaged over periods from 2 to 12 hours, which is a key quantity for atmospheric coupling by waves in Fig. 18.This plot demonstrates, that for these measurements the vertical flux of horizontal momentum rapidly decreases with altitude up to ∼45 km.Above ∼42 km it remains rather constant up to ∼70 km.In the altitude range from 42 to 70 km where we find a low variability of the momentum flux we analyzed its dependence on the horizontal propagation direction of the waves.The result is shown as polar histograms in Fig. 19.We see that the momentum flux of downward propagating waves is lower than that of upward propagating GW.Fig. 19 also shows that waves propagating nearly perpendicular to the mean wind carry the smallest flux for both up-and downward propagating GW.
Note, that the direction of the momentum flux is not necessarily along the major axis of the ellipse.The angle between the directions of momentum flux and GW propagation was estimated by Eq. 13 from Gavrilov et al. (1996) and does not exceed 2.8 • , which is much lower than the width of the bins in the histograms.

Summary and conclusion
In this paper a detailed step-by-step description of a new algorithm for derivation of GW parameters with justification for every step is presented.Most of these steps if considered independently, are well known and validated in numerous experimental works.The advantage and novelty of this work is their combination and some justifications of their importance and how they affect GW-analysis results.
Thus e.g., very first action normally performed on the measured time series is background removal.Since most conventional techniques based on smoothing or averaging in time or altitude ultimately introduce artifacts, we justify that application of the 2D-FFT for background removal is most appropriate.Advantage of this method is that it simultaneously accounts for both variability in space and time.
Specific feature of our algorithm for GW analysis is that it is insensitive to the particular background removal scheme.
Therefore, to avoid any degree of arbitrariness, the background removal can be excluded from fluctuation analysis when applying further steps of the analysis technique described in the manuscript.
As a next step we proposed to apply a scaling function of the form exp(z/(ςH)), where H is scale height, z is altitude, and the constant ς can be derived by a linear fit to fluctuation profiles and should be in a range 1-10 (we derived ς = 2.15 for our data).This, to our knowledge, a new technique which is not explicitly described in the literature.Advantage of this approach is to suppress exponential growth of GW-amplitudes to allow for equally weighted detection of wave signatures within the entire altitude range.This e.g., is clearly seen in wavelet scalograms which would otherwise be predominantly sensitive to strongest amplitudes, hiding out waves at lower altitudes.
The most essential part of the proposed analysis technique consist of fitting of cosines-waves to simultaneously measured profiles of winds and temperature and subsequent hodograph analysis of these fitted waves.We emphasize that this fit must be applied to all three quantities, i.e. zonal and meridional wind and temperature (u, v, and T), simultaneously.This ensures that we deal with a real GW which leaves its signature in all these physical quantities that were measured simultaneously in the same volume.The main difficulty in application of the hodograph analysis to real measurements is to find the wavelengths and altitude regions where certain GW dominates all measured quantities (u, v, and T).Since very often the measured data represent a mixture of vast different GWs, it is generally very difficult to find them automatically in the frame of hodograph analysis.Therefore, such work was always accomplished manually, by applying visual check of data and analysis quality.So was also done in particular by Baumgarten et al. (2015).The novelty of our approach is that its robustness ultimately allows for automation of the hodograph analysis.Also, our algorithm resolves many more GWs than it can be inferred by manually applied hodograph technique.
All these advantages are especially important since modern advanced measurement techniques (e.g.our lidar system described in Sec. 4) are capable of doing long duration measurements that cover large altitude range ∼30 to 80 km.This huge amount of data requires a robust and stable automatic analysis technique which we developed and presented in this work.
One obvious advantage of the proposed algorithm is that it allows for simultaneous detection of any kind of waves presented in the measurements.This includes not only GWs, but also tides.Since new analysis algorithm allows to apply a simplest background removal techniques like subtraction of mean, the necessity of removal tidal components a priori, which cannot be done unambiguously, is eliminated.All the detected waves can be sorted out on a statistical bases after the observational data base is analyzed by using the proposed algorithm.
Another specific feature of our analysis technique is the extension to the linear wave theory introduced in Sec. 3, the wave packet envelop term exp(−(z −z 0 ) 2 /2σ 2 ) that accounts for limited presence of the GW-packet in observations.This, however, only works in spatial domain, i.e. vertically.At the current stage of development our analysis technique is not capable of detecting life-time of gravity waves in observational data set.This capability is currently under development as well as an additional robust algorithm to pick out wave packets in time domain automatically.
By applying this new methodology to real data obtained by lidar during about 60 hours of observations in January 2016 we found 4507 single hodographs.In general, 5 to 10 waves were detected from every vertical profile.This allowed identifying and analyzing quasi monochromatic waves in about ∼80 % of the observations.The measurements were performed while a jet at the stratopause (45-55 km) of more than 100 m/s was located above the lidar station.We found a strong decrease in vertical flux of horizontal momentum up to ∼42 km altitude.Due to the strong wind above ∼40 km, it is likely that waves break, get absorbed, and reflected below this altitude region.The new method allows studying waves separated for up-and downward propagation according to their group velocities.
where w is vertical wind fluctuations and ρ is the atmospheric density.From continuity equation we get w = −(k /m) • u l and the vertical momentum flux is transformed to (e.g., Réchou et al., 2014): Author contributions.IS developed the analysis technique algorithm and code and performed the calculations.GB designed experiments and conducted measurements; GB and IS analyzed the data; IS, GB, and FJL contributed to the final manuscript.

Figure 1 .
Figure 1.Schematics of the method.(a) Altitude profile of horizontal velocity fluctuations.Blue dashed line demonstrates an envelope.Colored area marks altitude range of one wavelength where wave amplitude is most significant ([z0 − λz/2, z0 + λz/2]) (b) Hodograph ellipse of IGW horizontal velocity variations taken from altitude range marked in plot (a).Dashed line shows major axis of ellipse, which is a propagation direction of the wave.Numbers around ellipse are altitudes.In this schematics clockwise rotation

Figure 5 .
Figure 5. Wavelet transform of zonal and meridional wind and temperature fluctuations at 02:07 UT 10 January 2016.

Figure 9 .
Figure 9. Left: Number of waves detected per 1.5 km altitude range bin.Blue (orange) bars mark upward (downward) propagating GW.Middle: Mean coverage by detected waves when taking the altitude extent of the waves into account.The green profile indicates whether any wave was found, whereas blue and orange lines are for up-and downward propagating waves, respectively.Right: Background mean wind and temperature.

Figure 10 .
Figure 10.Reconstructed temperature fluctuations of upward propagating GW (left pannel) and downward propagating GW (right pannel).Contour lines show the background zonal wind and the numbers on the contour lines are given in m/s.

Figure 11 .
Figure 11.Color coded bars show the intrinsic period of upward (left) and downward (right) propagating waves.The length of the bar is given by the extend of the waves.Contour lines show the background total wind, the numbers on the contour lines are given in m/s.

Figure 12 .
Figure 12.Histograms of different GW properties separated for up and downward propagating waves.From left to right: vertical wavelength;intrinsic period; horizontal wavelength; horizontal phase speed.Estimated from equations 1, A5, A9, A10a, respectively.

Figure 13 .
Figure 13.Fourier power spectra of measured temperature fluctuations (blue) and of the reconstructed GWs (orange).

Figure 14 .Figure 15 .Figure 16 .Figure 17 .Figure 18 .
Figure 14.Kinetic (left) and potential (right) energies of waves.Colored boxes show 2-D histograms (number of waves per 1.5 km altitude, and 0.15 log(energy) bin).Lines show mean values of whole distribution (pink), upward/downward propagating waves (red/orange dashed) and black dashed lines are energies estimated from the variance of the temperature (left) and wind (right) fluctuations throughout the measurement.