A novel method for the extraction of local gravity wave parameters from gridded three-dimensional data: description, validation, and application

For the local diagnosis of wave properties, we develop, validate, and apply a novel method which is based on the Hilbert transform. It is called Unified Wave Diagnostics (UWaDi). It provides the wave amplitude and threedimensional wave number at any grid point for gridded threedimensional data. UWaDi is validated for a synthetic test case comprising two different wave packets. In comparison with other methods, the performance of UWaDi is very good with respect to wave properties and their location. For a first practical application of UWaDi, a minor sudden stratospheric warming on 30 January 2016 is chosen. Specifying the diagnostics for hydrostatic inertia–gravity waves in analyses from the European Centre for Medium-Range Weather Forecasts, we detect the local occurrence of gravity waves throughout the middle atmosphere. The local wave characteristics are discussed in terms of vertical propagation using the diagnosed local amplitudes and wave numbers. We also note some hints on local inertia–gravity wave generation by the stratospheric jet from the detection of shallow slow waves in the vicinity of its exit region.


Introduction
The importance of gravity waves (GWs) for the dynamics of the Earth's atmosphere is without controversy.They influence dynamics from planetary scales to turbulent microscales and play an important role in the middle atmosphere circulation (Fritts and Alexander, 2003).GWs typically appear as packets localised in space and time.Hence, it is desirable to diagnose them locally as precisely as possible.Here, we want to introduce a new method called Unified Wave Diagnosis (UWaDi).It provides phase-independent local wave quantities like amplitude and wave number without any prior assumption.In the following, we develop, validate, and apply this novel method.The application concentrates on the analysis of GWs for locally varying background wind conditions in the winter 2015-2016.
In the past, several methods were developed to estimate wave properties like amplitudes and wave number vectors.All of them have to deal with the fact that the data sampling procedure influences the results.A common approach to obtaining the vertical wave numbers and GW frequency of high-pass-filtered wind fluctuations is a method involving the Stokes parameters (Vincent and Fritts, 1987).This method is based on polarisation relations and works for single-column measurements.It provides the wave properties in preselected vertical height sections of finite lengths.Next to its original application to radar measurements, it is used for radiosonde data (Kramer et al., 2015).A supplement to this method called DIV was introduced by (Zülicke and Peters, 2006).It determines the dominating harmonic wave in a box from the first zero crossing of the autocorrelation function.The maximal detectable wavelength is restricted by the box size.The analysed quantity is the horizontal divergence to get the ageostrophic flow without numerical filtering.A further technique is based on sinusoidal few wave fits (S3-D) (Lehmann et al., 2012).This method was created for the analysis of three-dimensional data from remote sensing observations (Ern et al., 2017;Krisch et al., 2017) but is also applicable to model data (Preusse et al., 2014).The first modes L. Schoon and C. Zülicke: Diagnosis of local gravity wave properties with the highest variance are taken from a fit that minimises the variance-weighted squared deviations over all points in a finite analysis box.Only a small number of sinusoidal curves are fitted and there might remain uncovered variances in the analysis volume.These methods have in common that the analysed spatial scales are dependent on the predefined analysis box size and the assumption of the spatial homogeneity of the wave field in these boxes is essential.Nevertheless, these methods are superior to a classic Fourier transform in that they allow us to search for waves with bigger wavelengths than the box size.
Another three-dimensional spectral analysis method is the 3D Stockwell transform (3-D ST) (Wright et al., 2017).This method is capable of analysing the full range of length scales sampled in satellite data and is not restricted to fixed box sizes.At every grid point, a local wave spectrum is estimated using a window function of frequency-dependent width.S3-D and 3-D ST assume homogeneity inside the analysis box.Furthermore, they use a small number of the most prominent waves for the estimation of variances.This leaves some variance unattributed and hence means a loss of information.We search for a method which detects the full variance in each data point.
With UWaDi we find the dominating wave with the Hilbert transform at every data point.This approach does not rely on choosing the size of an analysis volume aforehand.The calculation of wave quantities at every grid point is computationally cheap.There is no need to assume homogeneity and no restriction on detectable wavelengths besides the Nyquist wavelength.Here, the method is developed to work with three-dimensional equally gridded data.In general, the Hilbert transform can be applied to data of any dimensionality.Wave properties such as the amplitude and wave number are estimated phase independently, while all variance is attributed to one wave mode.Every variable including any kind of wave-like structure can be diagnosed.(Zimin et al., 2003) used the method to obtain the envelope of a train of Rossby waves in one dimension.A supplement was made for waves not in line with grids by an extension of the formulation to stream lines to obtain quasi-one-dimensional wave packets (Zimin et al., 2006).(Sato et al., 2013) provide a threedimensional application to Rossby waves and GWs.(Glatt and Wirth, 2014) use the Hilbert transform to identify Rossby wave trains on a longitude-time plane and introduce their approach as an "objective identification method".Our method focuses on the local site of GW occurrence and the additional provision of the wave number in every dimension.The latter was not presented before.We aim to cover the retrieval of local wave properties from arbitrarily orientated wave packets.Amplitude and wave number are returned on the same grid as the input data.After the mathematical description of the method and how it is implemented, it will be validated with synthetic data to demonstrate its quality in comparison with other methods.
For a demonstration of a practical application in a geophysical context, we will investigate GWs.Their sources are usually found in the troposphere where waves are generated by flow over orography, by convection, frontal systems, and jet imbalances.These waves propagate upwards with increasing amplitudes and break in the middle atmosphere where they deposit their momentum to the background flow.Strong influence is exerted on global circulation patterns in the mesosphere and in the stratosphere (Holton, 1983;Garcia and Solomon, 1985).GWs play crucial roles in modulating the quasi-biennial oscillation (QBO) and the Brewer-Dobson circulation (Dunkerton, 1997;Alexander and Vincent, 2000;Ern et al., 2014).Another stratospheric phenomenon in which GWs play a role is a sudden stratospheric warming (SSW).For this phenomenon, a variety of definitions exist (Butler et al., 2015), but the most common one is given by the World Meteorological Organization stating that an SSW is characterised by a reversal of the 60 to 90 • N temperature gradient.Major warmings are associated with a wind reversal at 10 hPa and 60 • N, and minor SSWs (mSSWs) are associated with a wind deceleration at 10 hPa and 60 • N, where the prevailing westerlies are not turned into easterlies.While planetary waves are the most important drivers of SSWs (Andrews et al., 1987), GWs are affected by the differing background wind conditions during SSWs and are suspected to modulate the polar vortex in the post-warming phase of an SSW (Albers and Birner, 2014).The behaviour of GWs and planetary waves during an SSW was investigated through simulations and different measurement techniques.It was found that these GWs are, next to selective transmission, the subject of variable sources including unbalanced flow adjustment (Yamashita et al., 2010;Limpasuvan et al., 2011).We are interested in the longitudedependent transmission of GWs during an SSW.Pioneering work was done by (Dunkerton and Butchart, 1984).They analysed model data and found that the selective transmission of GWs during an SSW is dependent on longitude according to the planetary wave structures.Therefore, regions where vertical wave propagation is inhibited exist, as do regions where waves can propagate up to the mesosphere.The analysis of (Dunkerton and Butchart, 1984) was restricted to parameterised GWs of the "intermediate range" that they defined between 50 and 200 km.They state that it remains unclear in what way GWs of larger scale will act during SSWs.A study on a self-generated SSW in a model showed that GWs reverse the circulation in the mesosphere and lower thermosphere during an SSW by altering the altitude of GW breaking.This altitude is highly dependent on the specification of GW momentum flux in the lower atmosphere (Liu and Roble, 2002;Zülicke and Becker, 2013).In the UWaDi application, we want to locate the occurrence of GWs precisely in space and give first interpretations using the information on their changing amplitude and wave number in local vertical profiles.
The northern winter 2015-2016 brought up several interesting features, including specific GW patterns.The beginning of the winter was characterised by an extraordinarily strong and cold polar vortex driven by a deceleration of planetary waves in November-December 2015 (Matthias et al., 2016).Thereinafter, for the end of that winter a record Arctic ozone loss was expected (Manney and Lawrence, 2016).Furthermore, the extraordinary polar vortex caused a southward shift of planetary waves leading to anomalies in the QBO (Coy et al., 2017).A joint field campaign of the research projects METROSI, GW-LCYCLE2 (both part of ROMIC), and PACOG (MS-GWaves) took place in Scandinavia in January 2016.(Stober et al., 2017) found a summer-like zonal wind reversal in the upper mesosphere lasting until the end of January 2016, leading to different GW filtering processes in the mesosphere compared to usual winter-like wind conditions.During the field campaign the first tomographic observations of GWs by an infrared limb imager provided a full three-dimensional picture of a GW packet above Iceland (Krisch et al., 2017).Additionally, a remarkable comparative study indicates that forecasts of the current operational cycle (41r2) of the European Centre for Medium-Range Weather Forecasts (ECMWF) Integrated Forecast System (IFS) shows good accordance with space-borne lidar measurements while picturing large-scale and mesoscale wave structures in polar stratospheric clouds (Dörnbrack et al., 2017).We choose the mid-winter of 2016 for a first application of UWaDi because it is very well sampled with observations of GW properties.We intend to provide additional impulses to the evaluation of observations and start with a study of ECMWF analyses.These gridded data are suitable to analyse the local occurrence and their coupling as the analyses resolve essential parts of GW dynamics in the stratosphere.Validation studies with satellite measurements point out that ECMWF analyses capture GWs well in the middle and high latitudes (Yamashita et al., 2010;Preusse et al., 2014).Mid-latitude GWs are captured especially well, being driven by orographic and jet-stream-associated sources (Shutts and Vosper, 2011;Jewtoukoff et al., 2015).Our approach concentrates on fields of horizontal divergence in ECMWF IFS data.By choosing horizontal divergence as the analysis variable, the calculation of residuals is omitted.Hence, an explicit separation of a slowly varying background wind, which may include further uncertainties, is avoided.The horizontal divergence is a dynamical indicator for GWs because it is free from geostrophically balanced structures (Plougonven et al., 2003;Zülicke and Peters, 2006).Studying mountain waves, (Dörnbrack et al., 2012) and (Khaykin et al., 2015) have shown that divergence values correlate with GW-induced temperature anomalies.Furthermore, we use horizontal divergence to derive the total wave energy in Appendix B from the polarisation relations.
In this study we concentrate on vertical profiles of GW occurrence to give a first impression of the functionality of this method.Recently, the importance of the oblique propagation of GWs in a general context was discussed by (Yamashita et al., 2013), (Kalisch et al., 2014), and (Ehard et al., 2017).We point out that the meridional propagation of GWs is important for the deposition of GW drag from the stratosphere up into the mesosphere.For instance, (Krisch et al., 2017) discuss the role of oblique propagation for the redistribution of selective transmission of GWs in the upper troposphere and lower stratosphere.These are phenomena which require a detailed analysis of the propagation of localised wave packets.
Here, the authors focus on the introduction of the novel method and give first preliminary scientific results from a demonstrative application.We show locally diagnosed GW properties and give some hints on physical interpretation.A full three-dimensional spatial-temporal analysis of GWs during the SSW 2015-2016 goes beyond the scope of this paper and will be made the subject of subsequent publications.
The paper is organised as follows.After providing a stepby-step introduction and validation of the novel method in Sect.2, we give a short overview of the estimation of wave quantities for synthetic data and describe the analysis data.In Sect. 3 we show our results for the application of the mSSW on 30 January 2016 for which we study longitude-dependent GW occurrence.The discussion of our results in Sect. 4 is followed by the "Summary and conclusion" section (Sect.5).

Method and data
In this section we develop and validate an algorithm to extract wave parameters from gridded three-dimensional data.For the local diagnosis of waves, phase-independent estimates of wave amplitudes and the wave vector are essential.For this, we employ the Hilbert transform (e.g.Von Storch and Zwiers, 2001).The Hilbert transform shifts any sinusoidal wave structure by a quarter phase, i.e. turning a sine into a cosine.By constructing a new complex-valued data series consisting of the original field as the real part and its Hilbert transform as the imaginary part, the absolute value is always the amplitude (square root of squared real and imaginary part).The amplitude is independent of the phase of the wave and the wavelength of the underlying oscillation and there is no need for any explicit fitting of a particular wave.In addition, the absolute wave number in all three dimensions is determined from the phase gradient.The only requirement for the method to work is that the data contain any harmonic components.Thus, it works universally for any given variable, and hence it was called Unified Wave Diagnostics.

Step-by-step outline of the method
In the following we introduce UWaDi in a step-by-step outline.Further, we validate it with a well-defined test wave packet in comparison with other methods.In general, UWaDi is a script package which allows the user to steer data prepro- where r θ = R cos(θ) is the latitude-dependent Earth radius (R = 6371 km, Earth radius; θ , latitude in • ) and γ denotes the resolution of the gridded input data (e.g.γ = 0.36 • ).
2. To retrieve vertical equidistant levels, the hybrid levels of the ECMWF data are firstly transformed to pressure levels.Secondly, these pressure levels are assigned to equidistant height levels.For this purpose, we assume hydrostatic conditions and consider the surface geopotential and pressure as well as temperature and humidity.Both steps are performed with the help of common functions provided in the NCAR command language (NCL).This might cause problems in areas of high orography and inside the planetary boundary layer.These areas are not considered in the following analysis.
3. The method can handle any kind of variable.For the present application, we choose horizontal divergence.While this quantity was available in the ECMWF analyses, other data sources might require its calculation from the wind fields.However, the required preprocessing of the target variable is done in this step.
4. The underlying Hilbert transform is implemented with a discrete Fourier transform (DFT), which creates a complex spectrum in wave number space from the real valued data in real space (e.g.Smith, 1997).The processing of the three-dimensional data in the (x, y, z) space is begun with the x direction.The mathematics behind the Hilbert transform is described briefly for a onedimensional function f x originating in real space: The index k denotes the wave number space, and x describes the function f in real space.
5. DFTs can be biased by variance leakage through side lobes in spectral space.Tapering methods abandon this but can smear out nearby wave numbers.A loss of absolute amplitude can be overcome by using normalised weights (e.g.Von Storch and Zwiers, 2001).For the present study, however, the best results were obtained by turning the taper off.
6.In wave number space a rectangular bandpass filter reduces the complex spectrum to the user-predefined wave number limits k min and k max .Here, we make sure that only waves of the considered range of wave numbers are used for the following analysis.
7. To get back from wave number space an inverse DFT is performed.
8. The constructed complex-valued function fx consists of the input data f x as the real part and the Hilbert transformed function H (f x ) as the imaginary part: It provides the amplitude a x (Schönwiese, 2013), and the phase estimate x , x = atan 9. The phase gradient is a measure of the wave number modulus: 10. Due to the finite character of the data series it may happen that high-frequency spurious fluctuations appear after the Hilbert transform.We damp them by applying a low-pass filter.We smooth over a number of grid points determined by the lower wave number limit k min .
11.The identification of outliers is taken care of by two different quality checks.Firstly, the amplitude and wave number are checked for at least a half-undamped wave.Therefore, the packet length l x is essential.It is calculated by correlation function R i : with i max = N −1 5 (Chatfield, 2016).This method goes back to (Zülicke and Peters, 2006).The quality check is then defined by the inequality Secondly, the retrieved signals are required to lie above the noise level of the input data.An empirical threshold c checks the amplitude for being valid considering the standard deviation of the input function σ (f x ): Empirically, we use c = 0.01.This idea follows (Glatt and Wirth, 2014).UWaDi uses a quality flag q = 1, which is set to false (q = 0) if at least one quality check is rejected.
12. Steps 4 to 11 are repeated for the other dimensions (y, z).
13. Amplitude and absolute wave number are saved on the same grid as the input data to create a full threedimensional analysis of local wave quantities.The amplitude is combined to a wave-number-weighted sum of the three spatial dimensions: The absolute wave number is determined by with d denoting the spatial index.
The method provides an exact measure of the amplitude in the sense of the sum of the squared amplitudes of the wave modes.The dominating wave number is the amplitudeweighted sum of all.Spectrally wide dynamics can cause a significant reduction of information (Appendix A).Applying UWaDi with several narrow bandpass limits would provide information on spectrally spread waves.The wave numbers are estimated as moduli, which means the three-dimensional wave number (±k x , ±k y , ±k z ) allows for eight possible directions.However, the method is recommended for the first guess of the dominant wave packet, including the derivation of the intrinsic frequency from the dispersion relation.

Validation of the method
For a comparison of wave characteristics obtained with different methods we choose the test case presented in (Zimin et al., 2003) (Fig. 1a).In this exercise, a couple of localised wave packets with the wave numbers 4 and 9 are given in one dimension on the interval [0, 4π ] by Here, the quality check (step 11) requires the amplitudes to exceed half of the sample standard deviation.
UWaDi based on a Hilbert transform is a continuous method working without any box parameter.3-D ST, S3-D, and DIV need box width parameters to be adapted to the corresponding scientific case.A compromise between accuracy in space and wave number has to be found.For DIV the box length is set to L DIV = 8.0, which covers the largest anticipated wavelength.For 1-D ST, which was modified from 3-D ST, two steering parameters have to be adapted to the present task.We find a box width factor of C 1DST L = 0.25 suiting our requirements.Next to that, 1-D ST provides a correction factor for the absolute amplitude value C 1DST

A
. In three-dimensional data analyses enormous amplitude reductions can occur.Artificially changing the amplitude impacts the variance conservation of the technique.This is why we choose C 1DST A = 1 for our example.For S1-D a fixed box length has to be determined in advance.We find L S1−D = 1.5 to give acceptable results.We note that an extension to three dimensions would add information and therefore accuracy.However, in order to realise comparability we stay with the strictly one-dimensional set-up.
The method showing the best agreement with the theoretical value is UWaDi (Fig. 1b).For the amplitude, both wave packets are clearly distinguishable and the maximum peaks are recovered exactly.1-D ST and S1-D rebuild the wave packet shape as well.The lack of an absolute amplitude value might be adjusted with empirical correction factors provided for 1-D ST.Nevertheless, the amplitudes of both wave packets differ from each other.DIV provides a smeared out wavepacket envelope.The wave number calculation is best for UWaDi (Fig. 1c).The high peaks at the beginning and end of the wave packets are sorted out by the quality check.DIV meets the right wave numbers as well, but does not cover the whole spatial range of the two wave packets.1-D ST and S1-D show small deviations from the expected values.Altogether, UWaDi shows the best agreement with the theoretical expectations.

Analysis data
ECMWF data from the IFS operational cycle 41r1 is chosen for this analysis.We performed comparison studies between IFS data on different grid sizes (0.1, 0.36, 1 • ).By considering our bandpass filter conditions we found reliable results for the 0.1 and 0.36 • grids.To compromise between computational costs and the stability of the results we decided that data with a resolution of 0.36 • (ca.40 km) meet our requirements.They are retrieved from a resolution of T511.We discuss resolved gravity waves of a horizontal scale between 100 and 1500 km.In the vertical direction we are interested in gravity waves within the wavelength limits of 1 to 15 km.These scales fulfil the assumption of hydrostatics and cover the range of mid-and low-frequency GWs (Guest et al., 2000).
Vertical-propagating GWs are damped in ECMWF IFS products from 10 hPa (≈ 30 km) upwards (ECMWF, 2016).At 10 hPa the stratospheric sponge starts and a damping of wave propagation is expected (Jablonowski and Williamson, 2011).The mesospheric sponge follows at 1 hPa, acting on the divergence and therefore directly on the GW properties.We restrict our analysis to a maximum altitude of 45 km and therefore follow the advice of (Yamashita et al., 2010).We interpolate model levels to equidistant height levels between 2 and 45 km with a distance of 500 m and provide initial scientific analysis for a snapshot on 30 January 2016 at 00:00 UTC, corresponding to a minor SSW.

Gravity-wave-specific quantities
From the diagnosed fields of amplitude and wave number we calculate the kinematic wave energy e and wave action A. In order to find the ageostrophic GW motion we analyse fields of horizontal divergence.The kinematic wave energy is derived from polarisation equations for GWs assuming hydrostatics (Zülicke and Peters, 2006) (Appendix B): In this formula we need information on the divergence variance and the horizontal wave number.Both are provided by UWaDi from the three-dimensional divergence field.
The wave action is a conserved quantity as long as the slowly varying wave packet does not interact with the mean flow, e.g. by dissipation, absorption, or breaking (Bretherton, 1966).Wave action is defined by putting the kinematic wave energy e in relation to the intrinsic (flow-relative) frequency ω: with ρ being the density.The intrinsic frequency ω is calculated with the dispersion relation in mid-and low- ρ e ω = constant, one can see the following.
-Density effect: e ∝ 1 ρ ∝ exp z H .The above derived energy undergoes an exponential increase according to the density with the scale height H in vertical direction z.
-Wind effect: e ∝ u h .This relation holds for a stationary horizontally homogeneous mean flow (u(z), v(z)), which implies the invariance of the horizontal wave number (k x = const and k y = const) along with the apparent (ground-based) frequency (ω The energy scaling is obtained with the invariance of the wave action for an upwind wave (α k → α u +π ) as e = A ρ ω → A ρ ω+k h u h due to the Doppler shift of the intrinsic frequency (Marks and Eckermann, 1995).
For the following analysis primarily wave action is used.

Results
A minor SSW occurred on 30 January 2016.Figure 2a shows the wind speed of the Northern Hemisphere at 10 hPa at 00:00 UTC.A vortex displacement from the pole is visible.The displaced vortex causes areas of strongly curved winds.The horizontal divergence as a measure of GWs shows high wave activity above two areas (Fig. 2b).Firstly, above northern Europe horizontal divergence is aligned cross-stream.Secondly, spiral-like patterns appear above eastern Siberia, corresponding to an area of a curved jet streak.UWaDi applied to the field of horizontal divergence provides GW amplitude and wave action (Fig. 2c, d).GW amplitudes show patterns in regions of strongly alternating horizontal diver-gence.The wave action shows the highest peak above northern Europe and lower values above eastern Siberia.
In zonal mean the horizontal wavelength varies between 120 and 200 km (Fig 3a).In the mid-stratosphere between 18 and 40 km of altitude the horizontal wavelength remains nearly constant.Thus, our assumption of a homogeneous background wind field is approximately valid in the midstratosphere.The vertical wavelength scales from 2.2 to 5.2 km.It increases throughout the whole atmospheric section with a slight change in gradient at the altitude of the tropopause (10 km).The decrease in vertical wavelength above 35 km of altitude is dubious.It occurs frequently, also for other temporal snapshots.We suspect an influence of the IFS sponge layer, but do not exclude an influence of the decreasing zonal wind, and therefore remain critical on interpretations at these altitudes.The horizontal phase speed in the wave direction ) is approximated for an upwind wave (c h → ĉh − u h ).It remains unchanged for waves propagating passively through a stationary horizontally homogeneous wind field.The zonal mean remains nearly constant with a value of 5 ms −1 below the tropopause and then steadily increases in the stratosphere (Fig. 3b).The indicated variations in horizontal wave number and phase speed contradict assumptions of a passive propagation through a stationary horizontally homogeneous wind field.
We next inspect local profiles in different background wind conditions.Longitude-height sections of zonal wind (Fig. 4a) and wave action (Fig. 4b) at 60 • N help to find the location of interesting vertical profiles.Three profiles are chosen that are representative for regions of similar filter conditions.The first profile (1) at 7.56 • E is chosen to be in a longitudinal range characterised by strong zonal eastward winds and lies in the deceleration area of the jet stream above northern Europe.Profile (2) is at 151.92 • E in the area of a strongly curved stratospheric jet associated with the displacement of the polar vortex.In Fig. 4a it is visible as a wind intrusion in the altitude range between 14 and 34 km.The wave action shows a peak in that height range (Fig. 4b).For comparison we take a third profile (3) at 240.12 • E in a region of low wind speed above Canada, which means weak tropospheric and weak stratospheric jets.
To highlight the advantage of a local wave analysis we show profiles at selected longitudinal positions (Fig. 5a).
During a local increase in wind speed above northern Europe the vertical profiles of (1) show that the zonal wind meanders around 50 m s −1 (Fig. 5b).In the stratosphere, the vertical wavelength is nearly constant with an average wavelength of 8 km, which is higher in the troposphere with 11.5 km (not shown).The wave action shows a high gradient changing from 10 4 to 10 3 kg m −1 s −1 and remains at this high level above 20 km.
Above eastern Siberia a displaced stratospheric jet streak appears jointly with high wave action (Fig. 4).The zonal wind vertical profile (2) shows this in a height range of 14 to 30 km with an increase from 5 to maximal 30 m s −1 (Fig. 5c).The wave action follows the structure of the zonal wind.Notably, peak GW activity takes place in the lower stratosphere, clearly above the tropospheric jet stream.At these altitudes, GWs with a vertical wavelength of 2 km are found and the horizontal wavelength is about 350 km (Fig. 6).The last set of vertical profiles is located in an area of low zonal winds (3) (Fig. 5d).In the troposphere eastward winds and in the middle stratosphere westward winds occur with magnitudes below 20 ms −1 .Above the altitude of the wind reversal the wave action remains constant.

Discussion
The topic of selective wave transmission was first modelled by (Dunkerton and Butchart, 1984).They highlighted the longitude-dependent gravity wave propagation during an SSW by focussing on the impact on the mesosphere.(Ern et al., 2016) further point out that selective filtering by anomalous winds during an SSW create a heavy impact on GW propagation through the whole atmosphere.They point out theoretically that during the upward propagation of GWs, these waves get attenuated or eliminated by distinct specifications of background flows.Here, we compare local vertical profiles of background wind and GW parameters from analysis data.
Comparing cases (1), (2), and (3) with respect to their wave action profiles we diagnose at 42 km values ranging over 3 orders of magnitude.This confirms the high spatial variability of GWs during SSWs.Also the shapes of the profiles differ clearly.One class is characterised by a steady decrease until 25 km and constance above (such as the zonal mean and cases 1 and 3).Another class of profiles shows a well-expressed peak in the stratosphere and a steady decrease above (case 2).The detailed analysis-related GW dynamics goes beyond the scope of this paper.However, some hypotheses are formulated.
In the high-wind case (1), showing the highest values of wave action and weak in the vertical wave number above northern Europe, we find the longest vertical wavelength of our study (8 km).These steep waves hint at an orographic GW caused by the eastward flow above the Scandinavian mountain ridge Kjølen.This is comparable to the findings of (Limpasuvan et al., 2011), who showed that during the SSW 2009 westward-propagating GW packets emanate from key topographical features around the polar edge and that these wave packets have long vertical wavelengths.Further detailed analyses of this GW packet are expected in upcoming publications according to the joint measurement campaign of ROMIC and MS-GWaves at, amongst others, Kiruna, Sweden (67 • N, 20 • E).
In the displaced stratospheric jet case (2) (Fig. 5c) we find a GW packet triggered off a curved stratospheric jet streak.Firstly explained by (Uccellini and Koch, 1987), jet exit re- gions in the troposphere are expected to emit GWs.The increase in wave action in the middle stratosphere according to the intrusion of westerlies seen in Fig. 4a and b leads to the assumption that the present feature is associated with the stratospheric jet.The horizontal divergence field supports this hypothesis with GW structures spiraling out of the curved jet above eastern Siberia.These features are comparable to the simulations of the troposphere by (Mirzaei et al., 2014), which resolved shallow near-inertia waves in jet exit regions.Further hints are the higher wave action and the smallest found wavelength (1.9 km) along the largest found horizontal wavelength (350 km).
In the low-wind case (3) the high wave action in the upper troposphere up to the height of the wind reversal at 23 km might be caused by orographically induced GWs due to the position in the lee of the Rocky Mountains.Above that, the overall lowest values of wave action are found, agreeing with measurements in that height range (Thurairajah et al., 2010).It is interesting to note that the shape of the wave action profile is similar to the zonal mean profile (compare Figs. 3a  and 5c), while the wind above 25 km goes into another direction.

Summary and conclusion
With UWaDi we provide a tool for the analysis of gridded three-dimensional data to estimate amplitude and wave number phase independently and locally.The method is based on a Hilbert transform with which the use of predefined analysis boxes is avoided.It returns an estimate for each data grid point but stays computationally cheap.With regard to the locality it clearly shows its advantages in a method comparison for a synthetic one-dimensional test case.Disadvantages may play a role when the wave spectrum is broad and the attribu-    tion of the variance to one dominant harmonic is not justified.The additional estimation of the wave numbers completes the elements of a wave packet description.There is an ambiguity in the sign of the wave number and in the direction of the wave vector, which is the case for all spatial analysis methods considered in this paper.Still, the method is recommended as a reliable local diagnosis of medium complexity.
For the analysis of gravity waves, we estimated wave energy and wave action from the horizontal divergence.This approach does not require an explicit numerical filtering, which is a practical advantage.Other methods for the analysis of unbalanced flow components are available, although more complicated (Mirzaei et al., 2017).While the chosen formula requires the variance (or squared amplitude) and wave numbers, the Hilbert transform method may also provide local estimates for more complex quantities as included in the combined Rossby wave and gravity wave diagnostics of (Kinoshita and Sato, 2013).For our study, which is fo- cused on hydrostatic inertia-GW occurrence, the specific approach provides reliable results at low complexity.
With the demonstrative analysis of the synoptic situation on 30 January 2016 we show the advantages of UWaDi: providing wave quantities on every grid point.Longitudedependent GW filter processes, known as selective wave transmission, can be diagnosed spatially in detail.Local vertical profiles show selective wave transmissions and generation processes.We found cases with a steady decrease in the wave action through the tropopause up to the midstratosphere and constant values above in contrast to a case with a strong peak in the lower stratosphere and a steady decrease above.The latter happened in an area where the wind field is affected by the mSSW, characterised by a curved jet stream exit region in the stratosphere; we discuss GW generation by spontaneous emission.The diagnosed long horizontal and short vertical wavelengths support this hypothesis.With the present method we plan to join the closer evaluation of observations and models with respect to local features of GW generation and propagation.

Figure 2 .
Figure 2. Synoptical situation of the Northern Hemisphere from ECMWF analysis at 10 hPa on 30 January 2016.Wind speed (a) and horizontal divergence (b).Gravity wave amplitude (c) and wave action (d).Circled numbers along 60 • N latitude indicate the positions of three vertical profiles for later analysis.

Figure 4 .
Figure 4. Zonal wind (a) and wave action (b) at 60 • N on 30 January 2016 in longitude-height section.Numbered vertical profiles for further analysis are highlighted.